{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/NeuromatchAcademy/course-content/blob/master/tutorials/W3D2_DynamicNetworks/W3D2_Tutorial2.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "# Neuromatch Academy: Week 2, Day 4, Tutorial 2\n",
    "# Neuronal Network Dynamics: Wilson-Cowan Model\n",
    "\n",
    "__Content creators:__ Qinglong Gu, Songtin Li, Arvind Kumar, John Murray, Julijana Gjorgjieva \n",
    "\n",
    "__Content reviewers:__ Maryam Vaziri-Pashkam, Ella Batty, Lorenzo Fontolan, Richard Gao, Spiros Chavlis, Michael Waskom\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Tutorial Objectives\n",
    "In the previous tutorial, you became familiar with a neuronal network consisting of only an excitatory population. Here, we extend the approach we used to include both excitatory and inhibitory neuronal populations in our network. A simple, yet powerful model to study the dynamics of two interacting populations of excitatory and inhibitory neurons, is the so-called **Wilson-Cowan** rate model, which will be the subject of this tutorial.\n",
    "\n",
    "The objectives of this tutorial are to:\n",
    "\n",
    "- Write the **Wilson-Cowan** equations for the firing rate dynamics of a 2D system composed of an excitatory (E) and an inhibitory (I) population of neurons\n",
    "- Simulate the dynamics of the system, i.e., Wilson-Cowan model.\n",
    "- Plot the frequency-current (F-I) curves for both populations (i.e., E and I).\n",
    "- Visualize and inspect the behavior of the system using **phase plane analysis**, **vector fields**, and **nullclines**.\n",
    "\n",
    "Bonus steps:\n",
    "\n",
    "- Find and plot the **fixed points** of the Wilson-Cowan model.\n",
    "- Investigate the stability of the Wilson-Cowan model by linearizing its dynamics and examining the **Jacobian matrix**.\n",
    "- Learn how the Wilson-Cowan model can reach an oscillatory state.\n",
    "\n",
    "Bonus steps (applications):\n",
    "- Visualize the behavior of an Inhibition-stabilized network.\n",
    "- Simulate working memory using the Wilson-Cowan model.\n",
    "\n",
    "\\\\\n",
    "Reference paper:\n",
    "\n",
    "_[Wilson H and Cowan J (1972) Excitatory and inhibitory interactions in localized populations of model neurons. Biophysical Journal 12](https://doi.org/10.1016/S0006-3495(72)86068-5)_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:14.839905Z",
     "iopub.status.busy": "2021-05-25T01:14:14.839342Z",
     "iopub.status.idle": "2021-05-25T01:14:15.249704Z",
     "shell.execute_reply": "2021-05-25T01:14:15.250173Z"
    }
   },
   "outputs": [],
   "source": [
    "# Imports\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.optimize as opt       # root-finding algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:15.252974Z",
     "iopub.status.busy": "2021-05-25T01:14:15.252422Z",
     "iopub.status.idle": "2021-05-25T01:14:15.337491Z",
     "shell.execute_reply": "2021-05-25T01:14:15.336793Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title Figure Settings\n",
    "import ipywidgets as widgets       # interactive display\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "plt.style.use(\"https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/nma.mplstyle\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:15.355799Z",
     "iopub.status.busy": "2021-05-25T01:14:15.345323Z",
     "iopub.status.idle": "2021-05-25T01:14:15.370989Z",
     "shell.execute_reply": "2021-05-25T01:14:15.370523Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title Helper functions\n",
    "\n",
    "\n",
    "def default_pars(**kwargs):\n",
    "  pars = {}\n",
    "\n",
    "  # Excitatory parameters\n",
    "  pars['tau_E'] = 1.     # Timescale of the E population [ms]\n",
    "  pars['a_E'] = 1.2      # Gain of the E population\n",
    "  pars['theta_E'] = 2.8  # Threshold of the E population\n",
    "\n",
    "  # Inhibitory parameters\n",
    "  pars['tau_I'] = 2.0    # Timescale of the I population [ms]\n",
    "  pars['a_I'] = 1.0      # Gain of the I population\n",
    "  pars['theta_I'] = 4.0  # Threshold of the I population\n",
    "\n",
    "  # Connection strength\n",
    "  pars['wEE'] = 9.   # E to E\n",
    "  pars['wEI'] = 4.   # I to E\n",
    "  pars['wIE'] = 13.  # E to I\n",
    "  pars['wII'] = 11.  # I to I\n",
    "\n",
    "  # External input\n",
    "  pars['I_ext_E'] = 0.\n",
    "  pars['I_ext_I'] = 0.\n",
    "\n",
    "  # simulation parameters\n",
    "  pars['T'] = 50.        # Total duration of simulation [ms]\n",
    "  pars['dt'] = .1        # Simulation time step [ms]\n",
    "  pars['rE_init'] = 0.2  # Initial value of E\n",
    "  pars['rI_init'] = 0.2  # Initial value of I\n",
    "\n",
    "  # External parameters if any\n",
    "  for k in kwargs:\n",
    "      pars[k] = kwargs[k]\n",
    "\n",
    "  # Vector of discretized time points [ms]\n",
    "  pars['range_t'] = np.arange(0, pars['T'], pars['dt'])\n",
    "\n",
    "  return pars\n",
    "\n",
    "\n",
    "def F(x, a, theta):\n",
    "  \"\"\"\n",
    "  Population activation function, F-I curve\n",
    "\n",
    "  Args:\n",
    "    x     : the population input\n",
    "    a     : the gain of the function\n",
    "    theta : the threshold of the function\n",
    "\n",
    "  Returns:\n",
    "    f     : the population activation response f(x) for input x\n",
    "  \"\"\"\n",
    "\n",
    "  # add the expression of f = F(x)\n",
    "  f = (1 + np.exp(-a * (x - theta)))**-1 - (1 + np.exp(a * theta))**-1\n",
    "\n",
    "  return f\n",
    "\n",
    "\n",
    "def dF(x, a, theta):\n",
    "  \"\"\"\n",
    "  Derivative of the population activation function.\n",
    "\n",
    "  Args:\n",
    "    x     : the population input\n",
    "    a     : the gain of the function\n",
    "    theta : the threshold of the function\n",
    "\n",
    "  Returns:\n",
    "    dFdx  :  Derivative of the population activation function.\n",
    "  \"\"\"\n",
    "\n",
    "  dFdx = a * np.exp(-a * (x - theta)) * (1 + np.exp(-a * (x - theta)))**-2\n",
    "\n",
    "  return dFdx\n",
    "\n",
    "\n",
    "def plot_FI_inverse(x, a, theta):\n",
    "  f, ax = plt.subplots()\n",
    "  ax.plot(x, F_inv(x, a=a, theta=theta))\n",
    "  ax.set(xlabel=\"$x$\", ylabel=\"$F^{-1}(x)$\")\n",
    "\n",
    "\n",
    "def plot_FI_EI(x, FI_exc, FI_inh):\n",
    "  plt.figure()\n",
    "  plt.plot(x, FI_exc, 'b', label='E population')\n",
    "  plt.plot(x, FI_inh, 'r', label='I population')\n",
    "  plt.legend(loc='lower right')\n",
    "  plt.xlabel('x (a.u.)')\n",
    "  plt.ylabel('F(x)')\n",
    "  plt.show()\n",
    "\n",
    "\n",
    "def my_test_plot(t, rE1, rI1, rE2, rI2):\n",
    "\n",
    "  plt.figure()\n",
    "  ax1 = plt.subplot(211)\n",
    "  ax1.plot(pars['range_t'], rE1, 'b', label='E population')\n",
    "  ax1.plot(pars['range_t'], rI1, 'r', label='I population')\n",
    "  ax1.set_ylabel('Activity')\n",
    "  ax1.legend(loc='best')\n",
    "\n",
    "  ax2 = plt.subplot(212, sharex=ax1, sharey=ax1)\n",
    "  ax2.plot(pars['range_t'], rE2, 'b', label='E population')\n",
    "  ax2.plot(pars['range_t'], rI2, 'r', label='I population')\n",
    "  ax2.set_xlabel('t (ms)')\n",
    "  ax2.set_ylabel('Activity')\n",
    "  ax2.legend(loc='best')\n",
    "\n",
    "  plt.tight_layout()\n",
    "  plt.show()\n",
    "\n",
    "\n",
    "def plot_nullclines(Exc_null_rE, Exc_null_rI, Inh_null_rE, Inh_null_rI):\n",
    "\n",
    "  plt.figure()\n",
    "  plt.plot(Exc_null_rE, Exc_null_rI, 'b', label='E nullcline')\n",
    "  plt.plot(Inh_null_rE, Inh_null_rI, 'r', label='I nullcline')\n",
    "  plt.xlabel(r'$r_E$')\n",
    "  plt.ylabel(r'$r_I$')\n",
    "  plt.legend(loc='best')\n",
    "  plt.show()\n",
    "\n",
    "\n",
    "def my_plot_nullcline(pars):\n",
    "  Exc_null_rE = np.linspace(-0.01, 0.96, 100)\n",
    "  Exc_null_rI = get_E_nullcline(Exc_null_rE, **pars)\n",
    "  Inh_null_rI = np.linspace(-.01, 0.8, 100)\n",
    "  Inh_null_rE = get_I_nullcline(Inh_null_rI, **pars)\n",
    "\n",
    "  plt.plot(Exc_null_rE, Exc_null_rI, 'b', label='E nullcline')\n",
    "  plt.plot(Inh_null_rE, Inh_null_rI, 'r', label='I nullcline')\n",
    "  plt.xlabel(r'$r_E$')\n",
    "  plt.ylabel(r'$r_I$')\n",
    "  plt.legend(loc='best')\n",
    "\n",
    "\n",
    "def my_plot_vector(pars, my_n_skip=2, myscale=5):\n",
    "  EI_grid = np.linspace(0., 1., 20)\n",
    "  rE, rI = np.meshgrid(EI_grid, EI_grid)\n",
    "  drEdt, drIdt = EIderivs(rE, rI, **pars)\n",
    "\n",
    "  n_skip = my_n_skip\n",
    "\n",
    "  plt.quiver(rE[::n_skip, ::n_skip], rI[::n_skip, ::n_skip],\n",
    "             drEdt[::n_skip, ::n_skip], drIdt[::n_skip, ::n_skip],\n",
    "             angles='xy', scale_units='xy', scale=myscale, facecolor='c')\n",
    "\n",
    "  plt.xlabel(r'$r_E$')\n",
    "  plt.ylabel(r'$r_I$')\n",
    "\n",
    "\n",
    "def my_plot_trajectory(pars, mycolor, x_init, mylabel):\n",
    "  pars = pars.copy()\n",
    "  pars['rE_init'], pars['rI_init'] = x_init[0], x_init[1]\n",
    "  rE_tj, rI_tj = simulate_wc(**pars)\n",
    "\n",
    "  plt.plot(rE_tj, rI_tj, color=mycolor, label=mylabel)\n",
    "  plt.plot(x_init[0], x_init[1], 'o', color=mycolor, ms=8)\n",
    "  plt.xlabel(r'$r_E$')\n",
    "  plt.ylabel(r'$r_I$')\n",
    "\n",
    "\n",
    "def my_plot_trajectories(pars, dx, n, mylabel):\n",
    "  \"\"\"\n",
    "  Solve for I along the E_grid from dE/dt = 0.\n",
    "\n",
    "  Expects:\n",
    "  pars    : Parameter dictionary\n",
    "  dx      : increment of initial values\n",
    "  n       : n*n trjectories\n",
    "  mylabel : label for legend\n",
    "\n",
    "  Returns:\n",
    "    figure of trajectory\n",
    "  \"\"\"\n",
    "  pars = pars.copy()\n",
    "  for ie in range(n):\n",
    "    for ii in range(n):\n",
    "      pars['rE_init'], pars['rI_init'] = dx * ie, dx * ii\n",
    "      rE_tj, rI_tj = simulate_wc(**pars)\n",
    "      if (ie == n-1) & (ii == n-1):\n",
    "          plt.plot(rE_tj, rI_tj, 'gray', alpha=0.8, label=mylabel)\n",
    "      else:\n",
    "          plt.plot(rE_tj, rI_tj, 'gray', alpha=0.8)\n",
    "\n",
    "  plt.xlabel(r'$r_E$')\n",
    "  plt.ylabel(r'$r_I$')\n",
    "\n",
    "\n",
    "def plot_complete_analysis(pars):\n",
    "  plt.figure(figsize=(7.7, 6.))\n",
    "\n",
    "  # plot example trajectories\n",
    "  my_plot_trajectories(pars, 0.2, 6,\n",
    "                       'Sample trajectories \\nfor different init. conditions')\n",
    "  my_plot_trajectory(pars, 'orange', [0.6, 0.8],\n",
    "                     'Sample trajectory for \\nlow activity')\n",
    "  my_plot_trajectory(pars, 'm', [0.6, 0.6],\n",
    "                     'Sample trajectory for \\nhigh activity')\n",
    "\n",
    "  # plot nullclines\n",
    "  my_plot_nullcline(pars)\n",
    "\n",
    "  # plot vector field\n",
    "  EI_grid = np.linspace(0., 1., 20)\n",
    "  rE, rI = np.meshgrid(EI_grid, EI_grid)\n",
    "  drEdt, drIdt = EIderivs(rE, rI, **pars)\n",
    "  n_skip = 2\n",
    "  plt.quiver(rE[::n_skip, ::n_skip], rI[::n_skip, ::n_skip],\n",
    "             drEdt[::n_skip, ::n_skip], drIdt[::n_skip, ::n_skip],\n",
    "             angles='xy', scale_units='xy', scale=5., facecolor='c')\n",
    "\n",
    "  plt.legend(loc=[1.02, 0.57], handlelength=1)\n",
    "  plt.show()\n",
    "\n",
    "\n",
    "def plot_fp(x_fp, position=(0.02, 0.1), rotation=0):\n",
    "  plt.plot(x_fp[0], x_fp[1], 'ko', ms=8)\n",
    "  plt.text(x_fp[0] + position[0], x_fp[1] + position[1],\n",
    "           f'Fixed Point1=\\n({x_fp[0]:.3f}, {x_fp[1]:.3f})',\n",
    "           horizontalalignment='center', verticalalignment='bottom',\n",
    "           rotation=rotation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "The helper functions included:\n",
    "\n",
    "- Parameter dictionary: `default_pars(**kwargs)`. You can use:\n",
    "  - `pars = default_pars()` to get all the parameters, and then you can execute `print(pars)` to check these parameters. \n",
    "  - `pars = default_pars(T=T_sim, dt=time_step)` to set a different simulation time and time step\n",
    "  - After `pars = default_pars()`, use `par['New_para'] = value` to add a new parameter with its value\n",
    "  - Pass to functions that accept individual parameters with `func(**pars)`\n",
    "- F-I curve: `F(x, a, theta)`\n",
    "- Derivative of the F-I curve: `dF(x, a, theta)`\n",
    "- Plotting utilities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Section 1: Wilson-Cowan model of excitatory and inhibitory populations\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:15.379028Z",
     "iopub.status.busy": "2021-05-25T01:14:15.378471Z",
     "iopub.status.idle": "2021-05-25T01:14:15.425571Z",
     "shell.execute_reply": "2021-05-25T01:14:15.425076Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title Video 1: Phase analysis of the Wilson-Cowan E-I model\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id=\"GCpQmh45crM\", width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtube.com/watch?v=\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Section 1.1: Mathematical description of the WC model\n",
    "\n",
    "Many of the rich dynamics recorded in the brain are generated by the interaction of excitatory and inhibitory subtype neurons. Here, similar to what we did in the previous tutorial, we will model two coupled populations of E and I neurons (**Wilson-Cowan** model). We can write two coupled differential equations, each representing the dynamics of the excitatory or inhibitory population:\n",
    "\n",
    "\\begin{align}\n",
    "\\tau_E \\frac{dr_E}{dt} &= -r_E + F_E(w_{EE}r_E -w_{EI}r_I + I^{\\text{ext}}_E;a_E,\\theta_E)\\\\\n",
    "\\tau_I \\frac{dr_I}{dt} &= -r_I + F_I(w_{IE}r_E -w_{II}r_I + I^{\\text{ext}}_I;a_I,\\theta_I)    \\qquad (1)\n",
    "\\end{align}\n",
    "\n",
    "$r_E(t)$ represents the average activation (or firing rate) of the excitatory population at time $t$, and $r_I(t)$ the activation (or firing rate) of the inhibitory population. The parameters $\\tau_E$ and $\\tau_I$ control the timescales of the dynamics of each population. Connection strengths are given by: $w_{EE}$ (E $\\rightarrow$ E), $w_{EI}$ (I $\\rightarrow$ E), $w_{IE}$ (E $\\rightarrow$ I), and $w_{II}$ (I $\\rightarrow$ I). The terms $w_{EI}$ and $w_{IE}$ represent connections from inhibitory to excitatory population and vice versa, respectively. The transfer functions (or F-I curves) $F_E(x;a_E,\\theta_E)$ and $F_I(x;a_I,\\theta_I)$ can be different for the excitatory and the inhibitory populations.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Exercise 1: Plot out the F-I curves for the E and I populations \n",
    "\n",
    "Let's first plot out the F-I curves for the E and I populations using the function defined above with default parameter values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:15.430532Z",
     "iopub.status.busy": "2021-05-25T01:14:15.429652Z",
     "iopub.status.idle": "2021-05-25T01:14:15.433822Z",
     "shell.execute_reply": "2021-05-25T01:14:15.434243Z"
    }
   },
   "outputs": [],
   "source": [
    "pars = default_pars()\n",
    "x = np.arange(0, 10, .1)\n",
    "\n",
    "print(pars['a_E'], pars['theta_E'])\n",
    "print(pars['a_I'], pars['theta_I'])\n",
    "\n",
    "###################################################################\n",
    "# TODO for students: compute and plot the F-I curve here          #\n",
    "# Note: aE, thetaE, aI and theta_I are in the dictionray 'pars'   #\n",
    "###################################################################\n",
    "\n",
    "# Compute the F-I curve of the excitatory population\n",
    "FI_exc = ...\n",
    "\n",
    "# Compute the F-I curve of the inhibitory population\n",
    "FI_inh = ...\n",
    "\n",
    "# Uncomment when you fill the (...)\n",
    "# plot_FI_EI(x, FI_exc, FI_inh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:15.440019Z",
     "iopub.status.busy": "2021-05-25T01:14:15.439440Z",
     "iopub.status.idle": "2021-05-25T01:14:15.702574Z",
     "shell.execute_reply": "2021-05-25T01:14:15.703008Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "pars = default_pars()\n",
    "x = np.arange(0, 10, .1)\n",
    "\n",
    "print(pars['a_E'], pars['theta_E'])\n",
    "print(pars['a_I'], pars['theta_I'])\n",
    "\n",
    "# Compute the F-I curve of the excitatory population\n",
    "FI_exc = F(x, pars['a_E'], pars['theta_E'])\n",
    "# Compute the F-I curve of the inhibitory population\n",
    "FI_inh = F(x, pars['a_I'], pars['theta_I'])\n",
    "\n",
    "# Uncomment when you fill the (...)\n",
    "with plt.xkcd():\n",
    "  plot_FI_EI(x, FI_exc, FI_inh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Section 1.2: Simulation scheme for the Wilson-Cowan model\n",
    "\n",
    "Equation $1$ can be integrated numerically. Using the Euler method, the dynamics of E and I populations can be simulated on a time-grid of stepsize $\\Delta t$. The updates for the activity of the excitatory and the inhibitory populations can be written as:\n",
    "\n",
    "\\begin{align}\n",
    "r_E[k+1] &= r_E[k] + \\Delta r_E[k]\\\\\n",
    "r_I[k+1] &= r_I[k] + \\Delta r_I[k] \n",
    "\\end{align}\n",
    "\n",
    "with the increments\n",
    "\n",
    "\\begin{align}\n",
    "\\Delta r_E[k] &= \\frac{\\Delta t}{\\tau_E}[-r_E[k] + F_E(w_{EE}r_E[k] -w_{EI}r_I[k] + I^{\\text{ext}}_E[k];a_E,\\theta_E)]\\\\\n",
    "\\Delta r_I[k] &= \\frac{\\Delta t}{\\tau_I}[-r_I[k] + F_I(w_{IE}r_E[k] -w_{II}r_I[k] + I^{\\text{ext}}_I[k];a_I,\\theta_I)] \n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Exercise 2: Numerically integrate the Wilson-Cowan equations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:15.711315Z",
     "iopub.status.busy": "2021-05-25T01:14:15.709970Z",
     "iopub.status.idle": "2021-05-25T01:14:15.711971Z",
     "shell.execute_reply": "2021-05-25T01:14:15.712389Z"
    }
   },
   "outputs": [],
   "source": [
    "def simulate_wc(tau_E, a_E, theta_E, tau_I, a_I, theta_I,\n",
    "                wEE, wEI, wIE, wII, I_ext_E, I_ext_I,\n",
    "                rE_init, rI_init, dt, range_t, **other_pars):\n",
    "  \"\"\"\n",
    "  Simulate the Wilson-Cowan equations\n",
    "\n",
    "  Args:\n",
    "    Parameters of the Wilson-Cowan model\n",
    "\n",
    "  Returns:\n",
    "    rE, rI (arrays) : Activity of excitatory and inhibitory populations\n",
    "  \"\"\"\n",
    "  # Initialize activity arrays\n",
    "  Lt = range_t.size\n",
    "  rE = np.append(rE_init, np.zeros(Lt - 1))\n",
    "  rI = np.append(rI_init, np.zeros(Lt - 1))\n",
    "  I_ext_E = I_ext_E * np.ones(Lt)\n",
    "  I_ext_I = I_ext_I * np.ones(Lt)\n",
    "\n",
    "  # Simulate the Wilson-Cowan equations\n",
    "  for k in range(Lt - 1):\n",
    "    ########################################################################\n",
    "    # TODO for students: compute drE and drI and remove the error\n",
    "    raise NotImplementedError(\"Student exercise: compute the change in E/I\")\n",
    "    ########################################################################\n",
    "    # Calculate the derivative of the E population\n",
    "    drE = ...\n",
    "    # Calculate the derivative of the I population\n",
    "    drI = ...\n",
    "\n",
    "    # Update using Euler's method\n",
    "    rE[k + 1] = rE[k] + drE\n",
    "    rI[k + 1] = rI[k] + drI\n",
    "\n",
    "  return rE, rI\n",
    "\n",
    "\n",
    "pars = default_pars()\n",
    "\n",
    "# Here are two trajectories with close intial values\n",
    "# Uncomment these lines to test your function\n",
    "# rE1, rI1 = simulate_wc(**default_pars(rE_init=.32, rI_init=.15))\n",
    "# rE2, rI2 = simulate_wc(**default_pars(rE_init=.33, rI_init=.15))\n",
    "# my_test_plot(pars['range_t'], rE1, rI1, rE2, rI2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:15.730542Z",
     "iopub.status.busy": "2021-05-25T01:14:15.725150Z",
     "iopub.status.idle": "2021-05-25T01:14:16.138794Z",
     "shell.execute_reply": "2021-05-25T01:14:16.138312Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "def simulate_wc(tau_E, a_E, theta_E, tau_I, a_I, theta_I,\n",
    "                wEE, wEI, wIE, wII, I_ext_E, I_ext_I,\n",
    "                rE_init, rI_init, dt, range_t, **other_pars):\n",
    "  \"\"\"\n",
    "  Simulate the Wilson-Cowan equations\n",
    "\n",
    "  Args:\n",
    "    Parameters of the Wilson-Cowan model\n",
    "\n",
    "  Returns:\n",
    "    rE, rI (arrays) : Activity of excitatory and inhibitory populations\n",
    "  \"\"\"\n",
    "  # Initialize activity arrays\n",
    "  Lt = range_t.size\n",
    "  rE = np.append(rE_init, np.zeros(Lt - 1))\n",
    "  rI = np.append(rI_init, np.zeros(Lt - 1))\n",
    "  I_ext_E = I_ext_E * np.ones(Lt)\n",
    "  I_ext_I = I_ext_I * np.ones(Lt)\n",
    "\n",
    "  # Simulate the Wilson-Cowan equations\n",
    "  for k in range(Lt - 1):\n",
    "\n",
    "    # Calculate the derivative of the E population\n",
    "    drE = dt / tau_E * (-rE[k] + F(wEE * rE[k] - wEI * rI[k] + I_ext_E[k],\n",
    "                                   a_E, theta_E))\n",
    "\n",
    "    # Calculate the derivative of the I population\n",
    "    drI = dt / tau_I * (-rI[k] + F(wIE * rE[k] - wII * rI[k] + I_ext_I[k],\n",
    "                                   a_I, theta_I))\n",
    "\n",
    "    # Update using Euler's method\n",
    "    rE[k + 1] = rE[k] + drE\n",
    "    rI[k + 1] = rI[k] + drI\n",
    "\n",
    "  return rE, rI\n",
    "\n",
    "\n",
    "pars = default_pars()\n",
    "\n",
    "# Here are two trajectories with close intial values\n",
    "rE1, rI1 = simulate_wc(**default_pars(rE_init=.32, rI_init=.15))\n",
    "rE2, rI2 = simulate_wc(**default_pars(rE_init=.33, rI_init=.15))\n",
    "\n",
    "with plt.xkcd():\n",
    "  my_test_plot(pars['range_t'], rE1, rI1, rE2, rI2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "The two plots above show the temporal evolution of excitatory ($r_E$, blue) and inhibitory ($r_I$, red) activity for two different sets of initial conditions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Interactive Demo: population trajectories with different initial values\n",
    "In this interactive demo, we will simulate the Wilson-Cowan model and plot the trajectories of each population. What happens to the E and I population trajectories with different initial conditions?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:16.162666Z",
     "iopub.status.busy": "2021-05-25T01:14:16.162108Z",
     "iopub.status.idle": "2021-05-25T01:14:16.409731Z",
     "shell.execute_reply": "2021-05-25T01:14:16.409233Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "\n",
    "# @markdown Make sure you execute this cell to enable the widget!\n",
    "\n",
    "\n",
    "def plot_EI_diffinitial(rE_init=0.0):\n",
    "\n",
    "  pars = default_pars(rE_init=rE_init, rI_init=.15)\n",
    "  rE, rI = simulate_wc(**pars)\n",
    "\n",
    "  plt.figure()\n",
    "  plt.plot(pars['range_t'], rE, 'b', label='E population')\n",
    "  plt.plot(pars['range_t'], rI, 'r', label='I population')\n",
    "  plt.xlabel('t (ms)')\n",
    "  plt.ylabel('Activity')\n",
    "  plt.legend(loc='best')\n",
    "  plt.show()\n",
    "\n",
    "\n",
    "_ = widgets.interact(plot_EI_diffinitial, rE_init=(0.30, 0.35, .01))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Think!\n",
    "It is evident that the steady states of the neuronal response can be different when different initial states are chosen. Why is that? \n",
    "\n",
    "We will discuss this in the next section but try to think about it first."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Section 2: Phase plane analysis\n",
    "Just like we used a graphical method to study the dynamics of a 1-D system in the previous tutorial, here we will learn a  graphical approach called **phase plane analysis** to study the dynamics of a 2-D system like the Wilson-Cowan model. \n",
    "\n",
    "So far, we have plotted the activities of the two populations as a function of time, i.e., in the `Activity-t` plane, either the $(t, r_E(t))$ plane or the $(t, r_I(t))$ one. Instead, we can plot the two activities $r_E(t)$ and $r_I(t)$ against each other at any time point $t$. This characterization in the `rI-rE` plane $(r_I(t), r_E(t))$ is called the **phase plane**. Each line in the phase plane indicates how both $r_E$ and $r_I$ evolve with time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:16.416350Z",
     "iopub.status.busy": "2021-05-25T01:14:16.415766Z",
     "iopub.status.idle": "2021-05-25T01:14:16.458646Z",
     "shell.execute_reply": "2021-05-25T01:14:16.459104Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title Video 2: Nullclines and Vector Fields\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id=\"V2SBAK2Xf8Y\", width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtube.com/watch?v=\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Interactive Demo: From the Activity - time plane to the **$r_I$ - $r_E$** phase plane\n",
    "\n",
    "In this demo, we will visualize the system dynamics using both the `Activity-time` and the `(rE, rI)` phase plane. The circles indicate the activities at a given time $t$, while the lines represent the evolution of the system for the entire duration of the simulation.\n",
    "\n",
    "Move the time slider to better understand how the top plot relates to the bottom plot. Does the bottom plot have explicit information about time? What information does it give us?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:16.524438Z",
     "iopub.status.busy": "2021-05-25T01:14:16.471131Z",
     "iopub.status.idle": "2021-05-25T01:14:17.036288Z",
     "shell.execute_reply": "2021-05-25T01:14:17.035840Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "\n",
    "# @markdown Make sure you execute this cell to enable the widget!\n",
    "\n",
    "pars = default_pars(T=10, rE_init=0.6, rI_init=0.8)\n",
    "rE, rI = simulate_wc(**pars)\n",
    "\n",
    "\n",
    "def plot_activity_phase(n_t):\n",
    "  plt.figure(figsize=(8, 5.5))\n",
    "  plt.subplot(211)\n",
    "  plt.plot(pars['range_t'], rE, 'b', label=r'$r_E$')\n",
    "  plt.plot(pars['range_t'], rI, 'r', label=r'$r_I$')\n",
    "  plt.plot(pars['range_t'][n_t], rE[n_t], 'bo')\n",
    "  plt.plot(pars['range_t'][n_t], rI[n_t], 'ro')\n",
    "  plt.axvline(pars['range_t'][n_t], 0, 1, color='k', ls='--')\n",
    "  plt.xlabel('t (ms)', fontsize=14)\n",
    "  plt.ylabel('Activity', fontsize=14)\n",
    "  plt.legend(loc='best', fontsize=14)\n",
    "\n",
    "  plt.subplot(212)\n",
    "  plt.plot(rE, rI, 'k')\n",
    "  plt.plot(rE[n_t], rI[n_t], 'ko')\n",
    "  plt.xlabel(r'$r_E$', fontsize=18, color='b')\n",
    "  plt.ylabel(r'$r_I$', fontsize=18, color='r')\n",
    "\n",
    "  plt.tight_layout()\n",
    "  plt.show()\n",
    "\n",
    "\n",
    "_ = widgets.interact(plot_activity_phase, n_t=(0, len(pars['range_t']) - 1, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:17.042164Z",
     "iopub.status.busy": "2021-05-25T01:14:17.041692Z",
     "iopub.status.idle": "2021-05-25T01:14:17.044419Z",
     "shell.execute_reply": "2021-05-25T01:14:17.044888Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove explanation\n",
    "\n",
    "\"\"\"\n",
    "Discussion:\n",
    "\n",
    "- Phase plane portraits allows us to visualize out of all possible states which\n",
    "states a system can take. For example in this demo you see that the among all\n",
    "possible pairs of values for r1 and r2, this system can take only a limited number\n",
    "of states (those that lie on the black line).\n",
    "- There are other things we can infer from the phase portraits e.g. fixed\n",
    "points, trajectory to the fixed points etc.\n",
    "- Explicit information about time is not visible in the phase portraits\n",
    "\"\"\";"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Section 2.1: Nullclines of the Wilson-Cowan Equations\n",
    "\n",
    "An important concept in the phase plane analysis is the \"nullcline\" which is defined as the set of points in the phase plane where the activity of one population (but not necessarily the other) does not change.\n",
    "\n",
    "In other words, the $E$ and $I$ nullclines of Equation $(1)$ are defined as the points where $\\displaystyle{\\frac{dr_E}{dt}}=0$, for the excitatory nullcline, or $\\displaystyle\\frac{dr_I}{dt}=0$ for the inhibitory nullcline. That is:\n",
    "\n",
    "\\begin{align}\n",
    "-r_E + F_E(w_{EE}r_E -w_{EI}r_I + I^{\\text{ext}}_E;a_E,\\theta_E) &= 0  \\qquad (2)\\\\[1mm]\n",
    "-r_I + F_I(w_{IE}r_E -w_{II}r_I + I^{\\text{ext}}_I;a_I,\\theta_I) &= 0    \\qquad (3)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Exercise 3: Compute the nullclines of the Wilson-Cowan model\n",
    "\n",
    "In the next exercise, we will compute and plot the nullclines of the E and I population.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "Along the nullcline of excitatory population Equation $2$, you can calculate the inhibitory activity by rewriting Equation $2$ into\n",
    "\n",
    "\\begin{align}\n",
    "r_I = \\frac{1}{w_{EI}}\\big{[}w_{EE}r_E - F_E^{-1}(r_E; a_E,\\theta_E) + I^{\\text{ext}}_E \\big{]}. \\qquad(4)\n",
    "\\end{align}\n",
    "\n",
    "where $F_E^{-1}(r_E; a_E,\\theta_E)$ is the inverse of the excitatory transfer function (defined below). Equation $4$ defines the $r_E$ nullcline."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "Along the nullcline of inhibitory population Equation $3$, you can calculate the excitatory activity by rewriting Equation $3$ into \n",
    "\\begin{align}\n",
    "r_E = \\frac{1}{w_{IE}} \\big{[} w_{II}r_I + F_I^{-1}(r_I;a_I,\\theta_I) - I^{\\text{ext}}_I \\big{]}.  \\qquad (5)  \n",
    "\\end{align}\n",
    "\n",
    "shere $F_I^{-1}(r_I; a_I,\\theta_I)$ is the inverse of the inhibitory transfer function (defined below). Equation $5$ defines the $I$ nullcline."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "Note that, when computing the nullclines with Equations 4-5, we also need to calculate the inverse of the transfer functions. \\\\\n",
    "\n",
    "The inverse of the sigmoid shaped **f-I** function that we have been using is:\n",
    "\n",
    "$$F^{-1}(x; a, \\theta) = -\\frac{1}{a} \\ln \\left[ \\frac{1}{x + \\displaystyle \\frac{1}{1+\\text{e}^{a\\theta}}} - 1 \\right] + \\theta \\qquad (6)$$\n",
    "\n",
    "The first step is to implement the inverse transfer function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:17.050604Z",
     "iopub.status.busy": "2021-05-25T01:14:17.049402Z",
     "iopub.status.idle": "2021-05-25T01:14:17.051251Z",
     "shell.execute_reply": "2021-05-25T01:14:17.051787Z"
    }
   },
   "outputs": [],
   "source": [
    "def F_inv(x, a, theta):\n",
    "  \"\"\"\n",
    "  Args:\n",
    "    x         : the population input\n",
    "    a         : the gain of the function\n",
    "    theta     : the threshold of the function\n",
    "\n",
    "  Returns:\n",
    "    F_inverse : value of the inverse function\n",
    "  \"\"\"\n",
    "\n",
    "  #########################################################################\n",
    "  # TODO for students: compute F_inverse\n",
    "  raise NotImplementedError(\"Student exercise: compute the inverse of F(x)\")\n",
    "  #########################################################################\n",
    "  # Calculate Finverse (ln(x) can be calculated as np.log(x))\n",
    "  F_inverse = ...\n",
    "\n",
    "  return F_inverse\n",
    "\n",
    "\n",
    "pars = default_pars()\n",
    "x = np.linspace(1e-6, 1, 100)\n",
    "\n",
    "# Uncomment the next line to test your function\n",
    "# plot_FI_inverse(x, a=1, theta=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:17.076352Z",
     "iopub.status.busy": "2021-05-25T01:14:17.071518Z",
     "iopub.status.idle": "2021-05-25T01:14:17.308636Z",
     "shell.execute_reply": "2021-05-25T01:14:17.309077Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "def F_inv(x, a, theta):\n",
    "  \"\"\"\n",
    "  Args:\n",
    "    x         : the population input\n",
    "    a         : the gain of the function\n",
    "    theta     : the threshold of the function\n",
    "\n",
    "  Returns:\n",
    "    F_inverse : value of the inverse function\n",
    "  \"\"\"\n",
    "\n",
    "  # Calculate Finverse (ln(x) can be calculated as np.log(x))\n",
    "  F_inverse = -1/a * np.log((x + (1 + np.exp(a * theta))**-1)**-1 - 1) + theta\n",
    "\n",
    "  return F_inverse\n",
    "\n",
    "pars = default_pars()\n",
    "x = np.linspace(1e-6, 1, 100)\n",
    "\n",
    "with plt.xkcd():\n",
    "  plot_FI_inverse(x, a=1, theta=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "Now you can compute the nullclines, using Equations 4-5:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:17.316426Z",
     "iopub.status.busy": "2021-05-25T01:14:17.315120Z",
     "iopub.status.idle": "2021-05-25T01:14:17.317061Z",
     "shell.execute_reply": "2021-05-25T01:14:17.317510Z"
    }
   },
   "outputs": [],
   "source": [
    "def get_E_nullcline(rE, a_E, theta_E, wEE, wEI, I_ext_E, **other_pars):\n",
    "  \"\"\"\n",
    "  Solve for rI along the rE from drE/dt = 0.\n",
    "\n",
    "  Args:\n",
    "    rE    : response of excitatory population\n",
    "    a_E, theta_E, wEE, wEI, I_ext_E : Wilson-Cowan excitatory parameters\n",
    "    Other parameters are ignored\n",
    "\n",
    "  Returns:\n",
    "    rI    : values of inhibitory population along the nullcline on the rE\n",
    "  \"\"\"\n",
    "  #########################################################################\n",
    "  # TODO for students: compute rI for rE nullcline and disable the error\n",
    "  raise NotImplementedError(\"Student exercise: compute the E nullcline\")\n",
    "  #########################################################################\n",
    "  # calculate rI for E nullclines on rI\n",
    "  rI = ...\n",
    "\n",
    "  return rI\n",
    "\n",
    "\n",
    "def get_I_nullcline(rI, a_I, theta_I, wIE, wII, I_ext_I, **other_pars):\n",
    "  \"\"\"\n",
    "  Solve for E along the rI from dI/dt = 0.\n",
    "\n",
    "  Args:\n",
    "    rI    : response of inhibitory population\n",
    "    a_I, theta_I, wIE, wII, I_ext_I : Wilson-Cowan inhibitory parameters\n",
    "    Other parameters are ignored\n",
    "\n",
    "  Returns:\n",
    "    rE    : values of the excitatory population along the nullcline on the rI\n",
    "  \"\"\"\n",
    "  #########################################################################\n",
    "  # TODO for students: compute rI for rE nullcline and disable the error\n",
    "  raise NotImplementedError(\"Student exercise: compute the I nullcline\")\n",
    "  #########################################################################\n",
    "  # calculate rE for I nullclines on rI\n",
    "  rE = ...\n",
    "\n",
    "  return rE\n",
    "\n",
    "\n",
    "pars = default_pars()\n",
    "Exc_null_rE = np.linspace(-0.01, 0.96, 100)\n",
    "Inh_null_rI = np.linspace(-.01, 0.8, 100)\n",
    "\n",
    "# Uncomment these lines to test your functions\n",
    "# Exc_null_rI = get_E_nullcline(Exc_null_rE, **pars)\n",
    "# Inh_null_rE = get_I_nullcline(Inh_null_rI, **pars)\n",
    "\n",
    "# plot_nullclines(Exc_null_rE, Exc_null_rI, Inh_null_rE, Inh_null_rI)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:17.344942Z",
     "iopub.status.busy": "2021-05-25T01:14:17.338074Z",
     "iopub.status.idle": "2021-05-25T01:14:17.586411Z",
     "shell.execute_reply": "2021-05-25T01:14:17.586824Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "def get_E_nullcline(rE, a_E, theta_E, wEE, wEI, I_ext_E, **other_pars):\n",
    "  \"\"\"\n",
    "  Solve for rI along the rE from drE/dt = 0.\n",
    "\n",
    "  Args:\n",
    "    rE    : response of excitatory population\n",
    "    a_E, theta_E, wEE, wEI, I_ext_E : Wilson-Cowan excitatory parameters\n",
    "    Other parameters are ignored\n",
    "\n",
    "  Returns:\n",
    "    rI    : values of inhibitory population along the nullcline on the rE\n",
    "  \"\"\"\n",
    "  # calculate rI for E nullclines on rI\n",
    "  rI = 1 / wEI * (wEE * rE - F_inv(rE, a_E, theta_E) + I_ext_E)\n",
    "\n",
    "  return rI\n",
    "\n",
    "\n",
    "def get_I_nullcline(rI, a_I, theta_I, wIE, wII, I_ext_I, **other_pars):\n",
    "  \"\"\"\n",
    "  Solve for E along the rI from dI/dt = 0.\n",
    "\n",
    "  Args:\n",
    "    rI    : response of inhibitory population\n",
    "    a_I, theta_I, wIE, wII, I_ext_I : Wilson-Cowan inhibitory parameters\n",
    "    Other parameters are ignored\n",
    "\n",
    "  Returns:\n",
    "    rE    : values of the excitatory population along the nullcline on the rI\n",
    "  \"\"\"\n",
    "  # calculate rE for I nullclines on rI\n",
    "  rE = 1 / wIE * (wII * rI + F_inv(rI, a_I, theta_I) - I_ext_I)\n",
    "\n",
    "  return rE\n",
    "\n",
    "\n",
    "pars = default_pars()\n",
    "Exc_null_rE = np.linspace(-0.01, 0.96, 100)\n",
    "Inh_null_rI = np.linspace(-.01, 0.8, 100)\n",
    "\n",
    "# Uncomment these lines to test your functions\n",
    "Exc_null_rI = get_E_nullcline(Exc_null_rE, **pars)\n",
    "Inh_null_rE = get_I_nullcline(Inh_null_rI, **pars)\n",
    "\n",
    "with plt.xkcd():\n",
    "  plot_nullclines(Exc_null_rE, Exc_null_rI, Inh_null_rE, Inh_null_rI)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "Note that by definition along the blue line in the phase-plane spanned by $r_E, r_I$, $\\displaystyle{\\frac{dr_E(t)}{dt}} = 0$, therefore, it is called a nullcline. \n",
    "\n",
    "That is, the blue nullcline divides the phase-plane spanned by $r_E, r_I$ into two regions: on one side of the nullcline $\\displaystyle{\\frac{dr_E(t)}{dt}} > 0$ and on the other side $\\displaystyle{\\frac{dr_E(t)}{dt}} < 0$.\n",
    "\n",
    "The same is true for the red line along which $\\displaystyle{\\frac{dr_I(t)}{dt}} = 0$. That is, the red nullcline divides the phase-plane spanned by $r_E, r_I$ into two regions: on one side of the nullcline $\\displaystyle{\\frac{dr_I(t)}{dt}} > 0$ and on the other side $\\displaystyle{\\frac{dr_I(t)}{dt}} < 0$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Section 2.2: Vector field\n",
    "\n",
    "How can the phase plane and the nullcline curves help us understand the behavior of the Wilson-Cowan model? \n",
    "\n",
    "The activities of the $E$ and $I$ populations $r_E(t)$ and $r_I(t)$ at each time point $t$ correspond to a single point in the phase plane, with coordinates $(r_E(t),r_I(t))$. Therefore, the time-dependent trajectory of the system can be described as a continuous curve in the phase plane, and the tangent vector to the trajectory, which is defined as the vector $\\bigg{(}\\displaystyle{\\frac{dr_E(t)}{dt},\\frac{dr_I(t)}{dt}}\\bigg{)}$, indicates the direction towards which the activity is evolving and how fast is the activity changing along each axis. In fact, for each point $(E,I)$ in the phase plane, we can compute the tangent vector $\\bigg{(}\\displaystyle{\\frac{dr_E}{dt},\\frac{dr_I}{dt}}\\bigg{)}$, which  indicates the behavior of the system when it traverses that point. \n",
    "\n",
    "The map of tangent vectors in the phase plane is called **vector field**. The behavior of any trajectory in the phase plane is determined by i) the initial conditions $(r_E(0),r_I(0))$, and ii) the vector field $\\bigg{(}\\displaystyle{\\frac{dr_E(t)}{dt},\\frac{dr_I(t)}{dt}}\\bigg{)}$.\n",
    "\n",
    "In general, the value of the vector field at a particular point in the phase plane is represented by an arrow. The orientation and the size of the arrow reflect the direction and the norm of the vector, respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Exercise 4: Compute and plot the vector field $\\displaystyle{\\Big{(}\\frac{dr_E}{dt}, \\frac{dr_I}{dt} \\Big{)}}$\n",
    "\n",
    "Note that\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{dr_E}{dt} &= [-r_E + F_E(w_{EE}r_E -w_{EI}r_I + I^{\\text{ext}}_E;a_E,\\theta_E)]\\frac{1}{\\tau_E}\\\\\n",
    "\\frac{dr_I}{dt} &= [-r_I + F_I(w_{IE}r_E -w_{II}r_I + I^{\\text{ext}}_I;a_I,\\theta_I)]\\frac{1}{\\tau_I}    \n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:17.592898Z",
     "iopub.status.busy": "2021-05-25T01:14:17.592355Z",
     "iopub.status.idle": "2021-05-25T01:14:17.595789Z",
     "shell.execute_reply": "2021-05-25T01:14:17.595338Z"
    }
   },
   "outputs": [],
   "source": [
    "def EIderivs(rE, rI,\n",
    "             tau_E, a_E, theta_E, wEE, wEI, I_ext_E,\n",
    "             tau_I, a_I, theta_I, wIE, wII, I_ext_I,\n",
    "             **other_pars):\n",
    "  \"\"\"Time derivatives for E/I variables (dE/dt, dI/dt).\"\"\"\n",
    "  ######################################################################\n",
    "  # TODO for students: compute drEdt and drIdt and disable the error\n",
    "  raise NotImplementedError(\"Student exercise: compute the vector field\")\n",
    "  ######################################################################\n",
    "  # Compute the derivative of rE\n",
    "  drEdt = ...\n",
    "  # Compute the derivative of rI\n",
    "  drIdt = ...\n",
    "\n",
    "  return drEdt, drIdt\n",
    "\n",
    "\n",
    "# Uncomment below to test your function\n",
    "# plot_complete_analysis(default_pars())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:17.619764Z",
     "iopub.status.busy": "2021-05-25T01:14:17.606796Z",
     "iopub.status.idle": "2021-05-25T01:14:18.220156Z",
     "shell.execute_reply": "2021-05-25T01:14:18.220631Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "def EIderivs(rE, rI,\n",
    "             tau_E, a_E, theta_E, wEE, wEI, I_ext_E,\n",
    "             tau_I, a_I, theta_I, wIE, wII, I_ext_I,\n",
    "             **other_pars):\n",
    "  \"\"\"Time derivatives for E/I variables (dE/dt, dI/dt).\"\"\"\n",
    "  # Compute the derivative of rE\n",
    "  drEdt = (-rE + F(wEE * rE - wEI * rI + I_ext_E, a_E, theta_E)) / tau_E\n",
    "  # Compute the derivative of rI\n",
    "  drIdt = (-rI + F(wIE * rE - wII * rI + I_ext_I, a_I, theta_I)) / tau_I\n",
    "\n",
    "  return drEdt, drIdt\n",
    "\n",
    "\n",
    "with plt.xkcd():\n",
    "  plot_complete_analysis(default_pars())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "\n",
    "The last phase plane plot shows us that: \n",
    "- Trajectories seem to follow the direction of the vector field\n",
    "- Different trajectories eventually always reach one of two points depending on the initial conditions. \n",
    "- The two points where the trajectories converge are the intersection of the two nullcline curves. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Think! \n",
    "\n",
    "There are, in total, three intersection points, meaning that the system has three fixed points.\n",
    "\n",
    "- One of the fixed points (the one in the middle) is never the final state of a trajectory. Why is that? \n",
    "- Why the arrows tend to get smaller as they approach the fixed points?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:18.225028Z",
     "iopub.status.busy": "2021-05-25T01:14:18.224470Z",
     "iopub.status.idle": "2021-05-25T01:14:18.229740Z",
     "shell.execute_reply": "2021-05-25T01:14:18.230147Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove explanation\n",
    "\n",
    "\"\"\"\n",
    "Discussion:\n",
    "\n",
    "1. Because the middle fixed point is unstable. Trajectories only point to stable\n",
    "fixed points.\n",
    "\n",
    "2. The slope of dr/dt determines the speed at which the system states will evolve.\n",
    "At the nullclines dr_e/dt = 0 and/or dr_i/dt = 0. That means as we move close to\n",
    "the nullclines the dr/dt becomes smaller and smaller. Therefore the system state\n",
    "evolves slower and slower as we approach the nullcline.\n",
    "\"\"\";"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Summary\n",
    "\n",
    "Congratulations! You have finished the second day of the last week of the neuromatch academy! Here, you learned how to simulate a rate based model consisting of excitatory and inhibitory population of neurons.\n",
    "\n",
    "In the last tutorial on dynamical neuronal networks you learned to:\n",
    "- Implement and simulate a 2D system composed of an E and an I population of neurons using the **Wilson-Cowan** model\n",
    "- Plot the frequency-current (F-I) curves for both populations\n",
    "- Examine the behavior of the system using phase **plane analysis**, **vector fields**, and **nullclines**.\n",
    "\n",
    "Do you have more time? Have you finished early? We have more fun material for you!\n",
    "\n",
    "Below are some, more advanced concepts on dynamical systems:\n",
    "\n",
    "- You will learn how to find the fixed points on such a system, and to investigate its stability by linearizing its dynamics and examining the **Jacobian matrix**.\n",
    "- You will see identify conditions under which the Wilson-Cowan model can exhibit oscillations.\n",
    "\n",
    "If you need even more, there are two applications of the Wilson-Cowan model:\n",
    "\n",
    "- Visualization of an Inhibition-stabilized network\n",
    "- Simulation of working memory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Bonus 1: Fixed points, stability analysis, and limit cycles in the Wilson-Cowan model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:18.241818Z",
     "iopub.status.busy": "2021-05-25T01:14:18.241268Z",
     "iopub.status.idle": "2021-05-25T01:14:18.285509Z",
     "shell.execute_reply": "2021-05-25T01:14:18.285055Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title Video 3: Fixed points and their stability\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id=\"jIx26iQ69ps\", width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtube.com/watch?v=\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Fixed Points of the E/I system\n",
    "\n",
    "Clearly, the intersection points of the two nullcline curves are the fixed points of the Wilson-Cowan model in Equation $(1)$. \n",
    "\n",
    "In the next exercise, we will find the coordinate of all fixed points for a given set of parameters.\n",
    "\n",
    "We'll make use of two functions, similar to ones we saw in the previous tutorial, which use a root-finding algorithm to find the fixed points of the system with Excitatory and Inhibitory populations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:18.293464Z",
     "iopub.status.busy": "2021-05-25T01:14:18.292923Z",
     "iopub.status.idle": "2021-05-25T01:14:18.295895Z",
     "shell.execute_reply": "2021-05-25T01:14:18.296326Z"
    }
   },
   "outputs": [],
   "source": [
    "# @markdown *Execute the cell to define `my_fp` and `check_fp`*\n",
    "\n",
    "def my_fp(pars, rE_init, rI_init):\n",
    "  \"\"\"\n",
    "  Use opt.root function to solve Equations (2)-(3) from initial values\n",
    "  \"\"\"\n",
    "\n",
    "  tau_E, a_E, theta_E = pars['tau_E'], pars['a_E'], pars['theta_E']\n",
    "  tau_I, a_I, theta_I = pars['tau_I'], pars['a_I'], pars['theta_I']\n",
    "  wEE, wEI = pars['wEE'], pars['wEI']\n",
    "  wIE, wII = pars['wIE'], pars['wII']\n",
    "  I_ext_E, I_ext_I = pars['I_ext_E'], pars['I_ext_I']\n",
    "\n",
    "  # define the right hand of wilson-cowan equations\n",
    "  def my_WCr(x):\n",
    "\n",
    "    rE, rI = x\n",
    "    drEdt = (-rE + F(wEE * rE - wEI * rI + I_ext_E, a_E, theta_E)) / tau_E\n",
    "    drIdt = (-rI + F(wIE * rE - wII * rI + I_ext_I, a_I, theta_I)) / tau_I\n",
    "    y = np.array([drEdt, drIdt])\n",
    "\n",
    "    return y\n",
    "\n",
    "  x0 = np.array([rE_init, rI_init])\n",
    "  x_fp = opt.root(my_WCr, x0).x\n",
    "\n",
    "  return x_fp\n",
    "\n",
    "\n",
    "def check_fp(pars, x_fp, mytol=1e-6):\n",
    "  \"\"\"\n",
    "  Verify (drE/dt)^2 + (drI/dt)^2< mytol\n",
    "\n",
    "  Args:\n",
    "    pars    : Parameter dictionary\n",
    "    fp      : value of fixed point\n",
    "    mytol   : tolerance, default as 10^{-6}\n",
    "\n",
    "  Returns :\n",
    "    Whether it is a correct fixed point: True/False\n",
    "  \"\"\"\n",
    "\n",
    "  drEdt, drIdt = EIderivs(x_fp[0], x_fp[1], **pars)\n",
    "\n",
    "  return drEdt**2 + drIdt**2 < mytol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Exercise 5: Find the fixed points of the Wilson-Cowan model\n",
    "\n",
    "From the above nullclines, we notice that the system features  three fixed points with the parameters we used. To find their coordinates, we need to choose proper initial value to give to the `opt.root` function inside of the function `my_fp` we just defined, since the algorithm can only find fixed points in the vicinity of the initial value. \n",
    "\n",
    "In this exercise, you will use the function `my_fp` to find each of the fixed points by varying the initial values. Note that you can choose the values near the intersections of the nullclines as the initial values to calculate the fixed points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:18.300412Z",
     "iopub.status.busy": "2021-05-25T01:14:18.299851Z",
     "iopub.status.idle": "2021-05-25T01:14:18.303055Z",
     "shell.execute_reply": "2021-05-25T01:14:18.303538Z"
    }
   },
   "outputs": [],
   "source": [
    "pars = default_pars()\n",
    "\n",
    "######################################################################\n",
    "# TODO: Provide initial values to calculate the fixed points\n",
    "# Check if x_fp's are the correct with the function check_fp(x_fp)\n",
    "# Hint: vary different initial values to find the correct fixed points                                #\n",
    "######################################################################\n",
    "\n",
    "# my_plot_nullcline(pars)\n",
    "\n",
    "# Find the first fixed point\n",
    "# x_fp_1 = my_fp(pars, ..., ...)\n",
    "# if check_fp(pars, x_fp_1):\n",
    "#   plot_fp(x_fp_1)\n",
    "\n",
    "# Find the second fixed point\n",
    "# x_fp_2 = my_fp(pars, ..., ...)\n",
    "# if check_fp(pars, x_fp_2):\n",
    "#   plot_fp(x_fp_2)\n",
    "\n",
    "# Find the third fixed point\n",
    "# x_fp_3 = my_fp(pars, ..., ...)\n",
    "# if check_fp(pars, x_fp_3):\n",
    "#   plot_fp(x_fp_3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:18.325595Z",
     "iopub.status.busy": "2021-05-25T01:14:18.322626Z",
     "iopub.status.idle": "2021-05-25T01:14:18.662439Z",
     "shell.execute_reply": "2021-05-25T01:14:18.662909Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "pars = default_pars()\n",
    "\n",
    "with plt.xkcd():\n",
    "  my_plot_nullcline(pars)\n",
    "\n",
    "  # Find the first fixed point\n",
    "  x_fp_1 = my_fp(pars, 0.1, 0.1)\n",
    "  if check_fp(pars, x_fp_1):\n",
    "    plot_fp(x_fp_1)\n",
    "\n",
    "  # Find the second fixed point\n",
    "  x_fp_2 = my_fp(pars, 0.3, 0.3)\n",
    "  if check_fp(pars, x_fp_2):\n",
    "    plot_fp(x_fp_2)\n",
    "\n",
    "  # Find the third fixed point\n",
    "  x_fp_3 = my_fp(pars, 0.8, 0.6)\n",
    "  if check_fp(pars, x_fp_3):\n",
    "    plot_fp(x_fp_3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Stability of a fixed point and eigenvalues of the Jacobian Matrix\n",
    "\n",
    "First, let's first rewrite the system $1$ as:\n",
    "\n",
    "\\begin{align}\n",
    "&\\frac{dr_E}{dt} = G_E(r_E,r_I)\\\\[0.5mm]\n",
    "&\\frac{dr_I}{dt} = G_I(r_E,r_I)\n",
    "\\end{align}\n",
    "where\n",
    "\n",
    "\\begin{align}\n",
    "&G_E(r_E,r_I) = \\frac{1}{\\tau_E} [-r_E + F_E(w_{EE}r_E -w_{EI}r_I + I^{\\text{ext}}_E;a,\\theta)]\\\\[1mm]\n",
    "&G_I(r_E,r_I) = \\frac{1}{\\tau_I} [-r_I + F_I(w_{IE}r_E -w_{II}r_I + I^{\\text{ext}}_I;a,\\theta)]\n",
    "\\end{align}\n",
    "\n",
    "By definition, $\\displaystyle\\frac{dr_E}{dt}=0$ and $\\displaystyle\\frac{dr_I}{dt}=0$ at each fixed point. Therefore, if the initial state is exactly at the fixed point, the state of the system will not change as time evolves. \n",
    "\n",
    "However, if the initial state deviates slightly from the fixed point, there are two possibilities\n",
    "the trajectory will be attracted back to the \n",
    "\n",
    "1.   The trajectory will be attracted back to the fixed point\n",
    "2.   The trajectory will diverge from the fixed point. \n",
    "\n",
    "These two possibilities define the type of fixed point, i.e., stable or unstable. Similar to the 1D system studied in the previous tutorial, the stability of a fixed point $(r_E^*, r_I^*)$ can be determined by linearizing the dynamics of the system (can you figure out how?). The linearization will yield a matrix of first-order derivatives called the Jacobian matrix:\n",
    "\n",
    " \\begin{equation}\n",
    "   J=\n",
    "  \\left[ {\\begin{array}{cc}\n",
    "   \\displaystyle{\\frac{\\partial}{\\partial r_E}}G_E(r_E^*, r_I^*) & \\displaystyle{\\frac{\\partial}{\\partial r_I}}G_E(r_E^*, r_I^*)\\\\[1mm]\n",
    "   \\displaystyle\\frac{\\partial}{\\partial r_E} G_I(r_E^*, r_I^*) & \\displaystyle\\frac{\\partial}{\\partial r_I}G_I(r_E^*, r_I^*) \\\\\n",
    "  \\end{array} } \\right] \\quad (7)\n",
    "\\end{equation}\n",
    "\n",
    "\\\\\n",
    "\n",
    "The eigenvalues of the Jacobian matrix calculated at the fixed point will determine whether it is a stable or unstable fixed point.\n",
    "\n",
    "\\\\\n",
    "\n",
    "We can now compute the derivatives needed to build the Jacobian matrix. Using the chain and product rules the derivatives for the excitatory population are given by:\n",
    "\n",
    "\\\\\n",
    "\n",
    "\\begin{align}\n",
    "&\\frac{\\partial}{\\partial r_E} G_E(r_E^*, r_I^*) = \\frac{1}{\\tau_E} [-1 + w_{EE} F_E'(w_{EE}r_E^* -w_{EI}r_I^* +  I^{\\text{ext}}_E;\\alpha_E, \\theta_E)] \\\\[1mm]\n",
    "&\\frac{\\partial}{\\partial r_I} G_E(r_E^*, r_I^*)= \\frac{1}{\\tau_E} [-w_{EI} F_E'(w_{EE}r_E^* -w_{EI}r_I^* +  I^{\\text{ext}}_E;\\alpha_E, \\theta_E)]\n",
    "\\end{align}\n",
    "\n",
    "\\\\\n",
    "\n",
    "The same applies to the inhibitory population. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Exercise 6: Compute the Jacobian Matrix for the Wilson-Cowan model\n",
    "\n",
    "Here, you can use `dF(x,a,theta)` defined in the `Helper functions` to calculate the derivative of the F-I curve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:18.669430Z",
     "iopub.status.busy": "2021-05-25T01:14:18.668853Z",
     "iopub.status.idle": "2021-05-25T01:14:18.672749Z",
     "shell.execute_reply": "2021-05-25T01:14:18.672233Z"
    }
   },
   "outputs": [],
   "source": [
    "def get_eig_Jacobian(fp,\n",
    "                     tau_E, a_E, theta_E, wEE, wEI, I_ext_E,\n",
    "                     tau_I, a_I, theta_I, wIE, wII, I_ext_I, **other_pars):\n",
    "  \"\"\"Compute eigenvalues of the Wilson-Cowan Jacobian matrix at fixed point.\"\"\"\n",
    "  # Initialization\n",
    "  rE, rI = fp\n",
    "  J = np.zeros((2, 2))\n",
    "\n",
    "  ###########################################################################\n",
    "  # TODO for students: compute J and disable the error\n",
    "  raise NotImplementedError(\"Student excercise: compute the Jacobian matrix\")\n",
    "  ###########################################################################\n",
    "  # Compute the four elements of the Jacobian matrix\n",
    "  J[0, 0] = ...\n",
    "  J[0, 1] = ...\n",
    "  J[1, 0] = ...\n",
    "  J[1, 1] = ...\n",
    "\n",
    "  # Compute and return the eigenvalues\n",
    "  evals = np.linalg.eig(J)[0]\n",
    "  return evals\n",
    "\n",
    "\n",
    "# Uncomment below to test your function when you get the correct fixed point\n",
    "# eig_1 = get_eig_Jacobian(x_fp_1, **pars)\n",
    "# eig_2 = get_eig_Jacobian(x_fp_2, **pars)\n",
    "# eig_3 = get_eig_Jacobian(x_fp_3, **pars)\n",
    "\n",
    "# print(eig_1, 'Stable point')\n",
    "# print(eig_2, 'Unstable point')\n",
    "# print(eig_3, 'Stable point')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:18.685633Z",
     "iopub.status.busy": "2021-05-25T01:14:18.684306Z",
     "iopub.status.idle": "2021-05-25T01:14:18.687856Z",
     "shell.execute_reply": "2021-05-25T01:14:18.688300Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "def get_eig_Jacobian(fp,\n",
    "                     tau_E, a_E, theta_E, wEE, wEI, I_ext_E,\n",
    "                     tau_I, a_I, theta_I, wIE, wII, I_ext_I, **other_pars):\n",
    "  \"\"\"Compute eigenvalues of the Wilson-Cowan Jacobian matrix at fixed point.\"\"\"\n",
    "  # Initialization\n",
    "  rE, rI = fp\n",
    "  J = np.zeros((2, 2))\n",
    "\n",
    "  # Compute the four elements of the Jacobian matrix\n",
    "  J[0, 0] = (-1 + wEE * dF(wEE * rE - wEI * rI + I_ext_E,\n",
    "                           a_E, theta_E)) / tau_E\n",
    "\n",
    "  J[0, 1] = (-wEI * dF(wEE * rE - wEI * rI + I_ext_E,\n",
    "                       a_E, theta_E)) / tau_E\n",
    "\n",
    "  J[1, 0] = (wIE * dF(wIE * rE - wII * rI + I_ext_I,\n",
    "                      a_I, theta_I)) / tau_I\n",
    "\n",
    "  J[1, 1] = (-1 - wII * dF(wIE * rE - wII * rI + I_ext_I,\n",
    "                           a_I, theta_I)) / tau_I\n",
    "\n",
    "  # Compute and return the eigenvalues\n",
    "  evals = np.linalg.eig(J)[0]\n",
    "  return evals\n",
    "\n",
    "\n",
    "eig_1 = get_eig_Jacobian(x_fp_1, **pars)\n",
    "eig_2 = get_eig_Jacobian(x_fp_2, **pars)\n",
    "eig_3 = get_eig_Jacobian(x_fp_3, **pars)\n",
    "\n",
    "print(eig_1, 'Stable point')\n",
    "print(eig_2, 'Unstable point')\n",
    "print(eig_3, 'Stable point')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "As is evident, the stable fixed points correspond to the negative eigenvalues, while unstable point corresponds to at least one positive eigenvalue."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "The sign of the eigenvalues is determined by the connectivity (interaction) between excitatory and inhibitory populations. \n",
    "\n",
    "Below we investigate the effect of $w_{EE}$ on the nullclines and the eigenvalues of the dynamical system. \n",
    "\n",
    "\\* _Critical change is referred to as **pitchfork bifurcation**_. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Effect of `wEE` on the nullclines and the eigenvalues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:18.735064Z",
     "iopub.status.busy": "2021-05-25T01:14:18.733173Z",
     "iopub.status.idle": "2021-05-25T01:14:18.928813Z",
     "shell.execute_reply": "2021-05-25T01:14:18.928357Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "\n",
    "# @markdown Make sure you execute this cell to see the plot!\n",
    "\n",
    "eig_1_M = []\n",
    "eig_2_M = []\n",
    "eig_3_M = []\n",
    "pars = default_pars()\n",
    "wEE_grid = np.linspace(6, 10, 40)\n",
    "my_thre = 7.9\n",
    "for wEE in wEE_grid:\n",
    "  x_fp_1 = [0., 0.]\n",
    "  x_fp_2 = [.4, .1]\n",
    "  x_fp_3 = [.8, .1]\n",
    "  pars['wEE'] = wEE\n",
    "  if wEE < my_thre:\n",
    "    x_fp_1 = my_fp(pars, x_fp_1[0], x_fp_1[1])\n",
    "    eig_1 = get_eig_Jacobian(x_fp_1, **pars)\n",
    "    eig_1_M.append(np.max(np.real(eig_1)))\n",
    "  else:\n",
    "    x_fp_1 = my_fp(pars, x_fp_1[0], x_fp_1[1])\n",
    "    eig_1 = get_eig_Jacobian(x_fp_1, **pars)\n",
    "    eig_1_M.append(np.max(np.real(eig_1)))\n",
    "\n",
    "    x_fp_2 = my_fp(pars, x_fp_2[0], x_fp_2[1])\n",
    "    eig_2 = get_eig_Jacobian(x_fp_2, **pars)\n",
    "    eig_2_M.append(np.max(np.real(eig_2)))\n",
    "\n",
    "    x_fp_3 = my_fp(pars, x_fp_3[0], x_fp_3[1])\n",
    "    eig_3 = get_eig_Jacobian(x_fp_3, **pars)\n",
    "    eig_3_M.append(np.max(np.real(eig_3)))\n",
    "\n",
    "eig_1_M = np.array(eig_1_M)\n",
    "eig_2_M = np.array(eig_2_M)\n",
    "eig_3_M = np.array(eig_3_M)\n",
    "\n",
    "plt.figure(figsize=(8, 5.5))\n",
    "plt.plot(wEE_grid, eig_1_M, 'ko', alpha=0.5)\n",
    "plt.plot(wEE_grid[wEE_grid >= my_thre], eig_2_M, 'bo', alpha=0.5)\n",
    "plt.plot(wEE_grid[wEE_grid >= my_thre], eig_3_M, 'ro', alpha=0.5)\n",
    "plt.xlabel(r'$w_{\\mathrm{EE}}$')\n",
    "plt.ylabel('maximum real part of eigenvalue')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Interactive Demo: Nullclines position in the phase plane changes with parameter values\n",
    "\n",
    "In this interactive widget, we will explore how the nullclines move for different values of the parameter $w_{EE}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:18.955782Z",
     "iopub.status.busy": "2021-05-25T01:14:18.945782Z",
     "iopub.status.idle": "2021-05-25T01:14:19.501002Z",
     "shell.execute_reply": "2021-05-25T01:14:19.500360Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "\n",
    "# @markdown Make sure you execute this cell to enable the widget!\n",
    "\n",
    "\n",
    "def plot_nullcline_diffwEE(wEE):\n",
    "  \"\"\"\n",
    "    plot nullclines for different values of wEE\n",
    "  \"\"\"\n",
    "\n",
    "  pars = default_pars(wEE=wEE)\n",
    "\n",
    "  # plot the E, I nullclines\n",
    "  Exc_null_rE = np.linspace(-0.01, .96, 100)\n",
    "  Exc_null_rI = get_E_nullcline(Exc_null_rE, **pars)\n",
    "\n",
    "  Inh_null_rI = np.linspace(-.01, .8, 100)\n",
    "  Inh_null_rE = get_I_nullcline(Inh_null_rI, **pars)\n",
    "\n",
    "  plt.figure(figsize=(12, 5.5))\n",
    "  plt.subplot(121)\n",
    "  plt.plot(Exc_null_rE, Exc_null_rI, 'b', label='E nullcline')\n",
    "  plt.plot(Inh_null_rE, Inh_null_rI, 'r', label='I nullcline')\n",
    "  plt.xlabel(r'$r_E$')\n",
    "  plt.ylabel(r'$r_I$')\n",
    "  plt.legend(loc='best')\n",
    "\n",
    "  plt.subplot(222)\n",
    "  pars['rE_init'], pars['rI_init'] = 0.2, 0.2\n",
    "  rE, rI = simulate_wc(**pars)\n",
    "  plt.plot(pars['range_t'], rE, 'b', label='E population', clip_on=False)\n",
    "  plt.plot(pars['range_t'], rI, 'r', label='I population', clip_on=False)\n",
    "  plt.ylabel('Activity')\n",
    "  plt.legend(loc='best')\n",
    "  plt.ylim(-0.05, 1.05)\n",
    "  plt.title('E/I activity\\nfor different initial conditions',\n",
    "            fontweight='bold')\n",
    "\n",
    "  plt.subplot(224)\n",
    "  pars['rE_init'], pars['rI_init'] = 0.4, 0.1\n",
    "  rE, rI = simulate_wc(**pars)\n",
    "  plt.plot(pars['range_t'], rE, 'b', label='E population', clip_on=False)\n",
    "  plt.plot(pars['range_t'], rI, 'r', label='I population', clip_on=False)\n",
    "  plt.xlabel('t (ms)')\n",
    "  plt.ylabel('Activity')\n",
    "  plt.legend(loc='best')\n",
    "  plt.ylim(-0.05, 1.05)\n",
    "\n",
    "  plt.tight_layout()\n",
    "  plt.show()\n",
    "\n",
    "\n",
    "_ = widgets.interact(plot_nullcline_diffwEE, wEE=(6., 10., .01))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:19.508243Z",
     "iopub.status.busy": "2021-05-25T01:14:19.507584Z",
     "iopub.status.idle": "2021-05-25T01:14:19.511169Z",
     "shell.execute_reply": "2021-05-25T01:14:19.510727Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove explanation\n",
    "\n",
    "\"\"\"\n",
    "Discussion:\n",
    "\n",
    "- For low values of wEE there is only one fixed point and it is stable so initial\n",
    "conditions do not matter and the system always converge to the only fixed point\n",
    "\n",
    "- For high values of wEE we have three fixed points of which two are stable and\n",
    "one is unstable (or saddle). Now it matter where the initial conditions are. If\n",
    "the initial conditions are in the attractor region os the high activity fixed\n",
    "point then the system will converge to that (the bottom example).\n",
    "\"\"\";"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "We can also investigate the effect of different $w_{EI}$, $w_{IE}$, $w_{II}$, $\\tau_{E}$, $\\tau_{I}$, and $I_{E}^{\\text{ext}}$ on the stability of fixed points. In addition, we can also consider the perturbation of the parameters of the gain curve $F(\\cdot)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Limit cycle - Oscillations\n",
    "\n",
    "For some values of interaction terms ($w_{EE}, w_{IE}, w_{EI}, w_{II}$ the eigenvalues can become complex. When at least one pair of eigenvalues is complex, oscillations arise. \n",
    "The stability of oscillations is determined by the real part of the eigenvalues (+ve real part oscillations will grow, -ve real part oscillations will die out). The size of the complex part determines the frequency of oscillations. \n",
    "\n",
    "For instance, if we use a different set of parameters, $w_{EE}=6.4$, $w_{EI}=4.8$, $w_{IE}=6.$, $w_{II}=1.2$, and $I_{E}^{\\text{ext}}=0.8$, then we shall observe that the E and I population activity start to oscillate! Please execute the cell below to check the oscillatory behavior. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:19.546357Z",
     "iopub.status.busy": "2021-05-25T01:14:19.521347Z",
     "iopub.status.idle": "2021-05-25T01:14:19.762465Z",
     "shell.execute_reply": "2021-05-25T01:14:19.763053Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "\n",
    "# @markdown Make sure you execute this cell to see the oscillations!\n",
    "\n",
    "pars = default_pars(T=100.)\n",
    "pars['wEE'], pars['wEI'] = 6.4, 4.8\n",
    "pars['wIE'], pars['wII'] = 6.0, 1.2\n",
    "pars['I_ext_E'] = 0.8\n",
    "pars['rE_init'], pars['rI_init'] = 0.25, 0.25\n",
    "\n",
    "rE, rI = simulate_wc(**pars)\n",
    "plt.figure(figsize=(8, 5.5))\n",
    "plt.plot(pars['range_t'], rE, 'b', label=r'$r_E$')\n",
    "plt.plot(pars['range_t'], rI, 'r', label=r'$r_I$')\n",
    "plt.xlabel('t (ms)')\n",
    "plt.ylabel('Activity')\n",
    "plt.legend(loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Exercise 7: Plot the phase plane\n",
    "\n",
    "We can also understand the oscillations of the population behavior using the phase plane. By plotting a set of trajectories with different initial states, we can see that these trajectories will move in a circle instead of converging to a fixed point. This circle is called \"limit cycle\" and shows the periodic oscillations of the $E$ and $I$ population behavior under some conditions.\n",
    "\n",
    "Try to plot the phase plane using the previously defined functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:19.804154Z",
     "iopub.status.busy": "2021-05-25T01:14:19.790414Z",
     "iopub.status.idle": "2021-05-25T01:14:20.116682Z",
     "shell.execute_reply": "2021-05-25T01:14:20.117505Z"
    }
   },
   "outputs": [],
   "source": [
    "pars = default_pars(T=100.)\n",
    "pars['wEE'], pars['wEI'] = 6.4, 4.8\n",
    "pars['wIE'], pars['wII'] = 6.0, 1.2\n",
    "pars['I_ext_E'] = 0.8\n",
    "\n",
    "plt.figure(figsize=(7, 5.5))\n",
    "my_plot_nullcline(pars)\n",
    "###############################################################################\n",
    "# TODO for students: plot phase plane: nullclines, trajectories, fixed point  #\n",
    "###############################################################################\n",
    "\n",
    "# Find the correct fixed point\n",
    "# x_fp_1 = my_fp(pars, ..., ...)\n",
    "# if check_fp(pars, x_fp_1):\n",
    "#   plot_fp(x_fp_1, position=(0, 0), rotation=40)\n",
    "\n",
    "\n",
    "my_plot_trajectories(pars, 0.2, 3,\n",
    "                     'Sample trajectories \\nwith different initial values')\n",
    "my_plot_vector(pars)\n",
    "\n",
    "plt.legend(loc=[1.01, 0.7])\n",
    "plt.xlim(-0.05, 1.01)\n",
    "plt.ylim(-0.05, 0.65)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:20.145260Z",
     "iopub.status.busy": "2021-05-25T01:14:20.137867Z",
     "iopub.status.idle": "2021-05-25T01:14:20.578276Z",
     "shell.execute_reply": "2021-05-25T01:14:20.578715Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "pars = default_pars(T=100.)\n",
    "pars['wEE'], pars['wEI'] = 6.4, 4.8\n",
    "pars['wIE'], pars['wII'] = 6.0, 1.2\n",
    "pars['I_ext_E'] = 0.8\n",
    "\n",
    "with plt.xkcd():\n",
    "  plt.figure(figsize=(7, 5.5))\n",
    "  my_plot_nullcline(pars)\n",
    "\n",
    "  # Find the correct fixed point\n",
    "  x_fp_1 = my_fp(pars, 0.8, 0.8)\n",
    "  if check_fp(pars, x_fp_1):\n",
    "    plot_fp(x_fp_1, position=(0, 0), rotation=40)\n",
    "\n",
    "  my_plot_trajectories(pars, 0.2, 3,\n",
    "                       'Sample trajectories \\nwith different initial values')\n",
    "\n",
    "  my_plot_vector(pars)\n",
    "\n",
    "  plt.legend(loc=[1.01, 0.7])\n",
    "  plt.xlim(-0.05, 1.01)\n",
    "  plt.ylim(-0.05, 0.65)\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Interactive Demo: Limit cycle and oscillations.\n",
    "\n",
    "From the above examples, the change of model parameters changes the shape of the nullclines and, accordingly, the behavior of the $E$ and $I$ populations from steady fixed points to oscillations. However, the shape of the nullclines is unable to fully determine the behavior of the network. The vector field also matters. To demonstrate this, here, we will investigate the effect of time constants on the population behavior. By changing the inhibitory time constant $\\tau_I$, the nullclines do not change, but the network behavior changes substantially from steady state to oscillations with different frequencies. \n",
    "\n",
    "Such a dramatic change in the system behavior is referred to as a **bifurcation**. \n",
    "\n",
    "\\\\\n",
    "Please execute the code below to check this out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:20.662437Z",
     "iopub.status.busy": "2021-05-25T01:14:20.661219Z",
     "iopub.status.idle": "2021-05-25T01:14:21.436088Z",
     "shell.execute_reply": "2021-05-25T01:14:21.435525Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "\n",
    "# @markdown Make sure you execute this cell to enable the widget!\n",
    "\n",
    "\n",
    "def time_constant_effect(tau_i=0.5):\n",
    "\n",
    "  pars = default_pars(T=100.)\n",
    "  pars['wEE'], pars['wEI'] = 6.4, 4.8\n",
    "  pars['wIE'], pars['wII'] = 6.0, 1.2\n",
    "  pars['I_ext_E'] = 0.8\n",
    "\n",
    "  pars['tau_I'] = tau_i\n",
    "\n",
    "  Exc_null_rE = np.linspace(0.0, .9, 100)\n",
    "  Inh_null_rI = np.linspace(0.0, .6, 100)\n",
    "\n",
    "  Exc_null_rI = get_E_nullcline(Exc_null_rE, **pars)\n",
    "  Inh_null_rE = get_I_nullcline(Inh_null_rI, **pars)\n",
    "\n",
    "  plt.figure(figsize=(12.5, 5.5))\n",
    "\n",
    "  plt.subplot(121)  # nullclines\n",
    "  plt.plot(Exc_null_rE, Exc_null_rI, 'b', label='E nullcline', zorder=2)\n",
    "  plt.plot(Inh_null_rE, Inh_null_rI, 'r', label='I nullcline', zorder=2)\n",
    "  plt.xlabel(r'$r_E$')\n",
    "  plt.ylabel(r'$r_I$')\n",
    "\n",
    "  # fixed point\n",
    "  x_fp_1 = my_fp(pars, 0.5, 0.5)\n",
    "  plt.plot(x_fp_1[0], x_fp_1[1], 'ko', zorder=2)\n",
    "\n",
    "  eig_1 = get_eig_Jacobian(x_fp_1, **pars)\n",
    "\n",
    "  # trajectories\n",
    "  for ie in range(5):\n",
    "    for ii in range(5):\n",
    "      pars['rE_init'], pars['rI_init'] = 0.1 * ie, 0.1 * ii\n",
    "      rE_tj, rI_tj = simulate_wc(**pars)\n",
    "      plt.plot(rE_tj, rI_tj, 'k', alpha=0.3, zorder=1)\n",
    "\n",
    "  # vector field\n",
    "  EI_grid_E = np.linspace(0., 1.0, 20)\n",
    "  EI_grid_I = np.linspace(0., 0.6, 20)\n",
    "  rE, rI = np.meshgrid(EI_grid_E, EI_grid_I)\n",
    "  drEdt, drIdt = EIderivs(rE, rI, **pars)\n",
    "  n_skip = 2\n",
    "  plt.quiver(rE[::n_skip, ::n_skip], rI[::n_skip, ::n_skip],\n",
    "              drEdt[::n_skip, ::n_skip], drIdt[::n_skip, ::n_skip],\n",
    "              angles='xy', scale_units='xy', scale=10, facecolor='c')\n",
    "  plt.title(r'$\\tau_I=$'+'%.1f ms' % tau_i)\n",
    "\n",
    "  plt.subplot(122)  # sample E/I trajectories\n",
    "  pars['rE_init'], pars['rI_init'] = 0.25, 0.25\n",
    "  rE, rI = simulate_wc(**pars)\n",
    "  plt.plot(pars['range_t'], rE, 'b', label=r'$r_E$')\n",
    "  plt.plot(pars['range_t'], rI, 'r', label=r'$r_I$')\n",
    "  plt.xlabel('t (ms)')\n",
    "  plt.ylabel('Activity')\n",
    "  plt.title(r'$\\tau_I=$'+'%.1f ms' % tau_i)\n",
    "  plt.legend(loc='best')\n",
    "  plt.tight_layout()\n",
    "  plt.show()\n",
    "\n",
    "\n",
    "_ = widgets.interact(time_constant_effect, tau_i=(0.2, 3, .1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "Both $\\tau_E$ and $\\tau_I$ feature in the Jacobian of the two population network (eq 7). So here is seems that the by increasing $\\tau_I$ the eigenvalues corresponding to the stable fixed point are becoming complex.\n",
    "\n",
    "Intuitively, when $\\tau_I$ is smaller, inhibitory activity changes faster than excitatory activity. As inhibition exceeds above a certain value, high inhibition inhibits excitatory population but that in turns means that inhibitory population gets smaller input (from the exc. connection). So inhibition decreases rapidly. But this means that excitation recovers -- and so on ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Bonus 2: Inhibition-stabilized network (ISN)\n",
    "\n",
    "As described above, one can obtain the linear approximation around the fixed point as \n",
    "\n",
    " \\begin{equation}\n",
    "   \\frac{d}{dr}  \\vec{R}=\n",
    "  \\left[ {\\begin{array}{cc}\n",
    "   \\displaystyle{\\frac{\\partial G_E}{\\partial r_E}} & \\displaystyle{\\frac{\\partial G_E}{\\partial r_I}}\\\\[1mm]\n",
    "   \\displaystyle\\frac{\\partial G_I}{\\partial r_E} & \\displaystyle\\frac{\\partial G_I}{\\partial r_I} \\\\\n",
    "  \\end{array} } \\right] \\vec{R},\n",
    "\\end{equation}\n",
    "\n",
    "\\\\\n",
    "\n",
    "\n",
    "where $\\vec{R} = [r_E, r_I]^{\\rm T}$ is the vector of the E/I activity.\n",
    "\n",
    "Let's direct our attention to the excitatory subpopulation which follows:\n",
    "\n",
    "\\\\\n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{dr_E}{dt} = \\frac{\\partial G_E}{\\partial r_E}\\cdot r_E + \\frac{\\partial G_E}{\\partial r_I} \\cdot r_I\n",
    "\\end{equation}\n",
    "\n",
    "\\\\\n",
    "\n",
    "Recall that, around fixed point $(r_E^*, r_I^*)$:\n",
    "\n",
    "\\\\\n",
    "\n",
    "\\begin{align}\n",
    "&\\frac{\\partial}{\\partial r_E}G_E(r_E^*, r_I^*) = \\frac{1}{\\tau_E} [-1 + w_{EE} F'_{E}(w_{EE}r_E^* -w_{EI}r_I^* + I^{\\text{ext}}_E; \\alpha_E, \\theta_E)] \\qquad (8)\\\\[1mm]\n",
    "&\\frac{\\partial}{\\partial r_I}G_E(r_E^*, r_I^*) = \\frac{1}{\\tau_E} [-w_{EI} F'_{E}(w_{EE}r_E^* -w_{EI}r_I^* + I^{\\text{ext}}_E; \\alpha_E, \\theta_E)] \\qquad (9)\\\\[1mm]\n",
    "&\\frac{\\partial}{\\partial r_E}G_I(r_E^*, r_I^*) = \\frac{1}{\\tau_I} [w_{IE} F'_{I}(w_{IE}r_E^* -w_{II}r_I^* + I^{\\text{ext}}_I; \\alpha_I, \\theta_I)] \\qquad (10)\\\\[1mm]\n",
    "&\\frac{\\partial}{\\partial r_I}G_I(r_E^*, r_I^*) = \\frac{1}{\\tau_I} [-1-w_{II} F'_{I}(w_{IE}r_E^* -w_{II}r_I^* + I^{\\text{ext}}_I; \\alpha_I, \\theta_I)] \\qquad (11)\n",
    "\\end{align} \\\\\n",
    "\n",
    "\n",
    "From Equation. (8), it is clear that $\\displaystyle{\\frac{\\partial G_E}{\\partial r_I}}$ is negative since the $\\displaystyle{\\frac{dF}{dx}}$ is always positive. It can be understood by that the recurrent inhibition from the inhibitory activity ($I$) can reduce the excitatory ($E$) activity. However, as described above, $\\displaystyle{\\frac{\\partial G_E}{\\partial r_E}}$ has negative terms related to the \"leak\" effect, and positive term related to the recurrent excitation. Therefore, it leads to two different regimes:\n",
    "\n",
    "- $\\displaystyle{\\frac{\\partial}{\\partial r_E}G_E(r_E^*, r_I^*)}<0$, **noninhibition-stabilized\n",
    "network (non-ISN) regime**\n",
    "\n",
    "- $\\displaystyle{\\frac{\\partial}{\\partial r_E}G_E(r_E^*, r_I^*)}>0$, **inhibition-stabilized\n",
    "network (ISN) regime**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Exercise 8: Compute $\\displaystyle{\\frac{\\partial G_E}{\\partial r_E}}$\n",
    "Implemet the function to calculate the $\\displaystyle{\\frac{\\partial G_E}{\\partial r_E}}$ for the default parameters, and the parameters of the limit cycle case. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:21.448010Z",
     "iopub.status.busy": "2021-05-25T01:14:21.446706Z",
     "iopub.status.idle": "2021-05-25T01:14:21.449821Z",
     "shell.execute_reply": "2021-05-25T01:14:21.449332Z"
    }
   },
   "outputs": [],
   "source": [
    "def get_dGdE(fp, tau_E, a_E, theta_E, wEE, wEI, I_ext_E, **other_pars):\n",
    "  \"\"\"\n",
    "  Simulate the Wilson-Cowan equations\n",
    "\n",
    "  Args:\n",
    "    fp   : fixed point (E, I), array\n",
    "    Other arguments are parameters of the Wilson-Cowan model\n",
    "\n",
    "  Returns:\n",
    "    J    : the 2x2 Jacobian matrix\n",
    "  \"\"\"\n",
    "  rE, rI = fp\n",
    "\n",
    "  ##########################################################################\n",
    "  # TODO for students: compute dGdrE and disable the error\n",
    "  raise NotImplementedError(\"Student excercise: compute the dG/dE, Eq. (13)\")\n",
    "  ##########################################################################\n",
    "  # Calculate the J[0,0]\n",
    "  dGdrE = ...\n",
    "\n",
    "  return dGdrE\n",
    "\n",
    "\n",
    "# Uncomment below to test your function\n",
    "pars = default_pars()\n",
    "x_fp_1 = my_fp(pars, 0.1, 0.1)\n",
    "x_fp_2 = my_fp(pars, 0.3, 0.3)\n",
    "x_fp_3 = my_fp(pars, 0.8, 0.6)\n",
    "\n",
    "# dGdrE1 = get_dGdE(x_fp_1, **pars)\n",
    "# dGdrE2 = get_dGdE(x_fp_2, **pars)\n",
    "# dGdrE3 = get_dGdE(x_fp_3, **pars)\n",
    "\n",
    "print(f'For the default case:')\n",
    "# print(f'dG/drE(fp1) = {dGdrE1:.3f}')\n",
    "# print(f'dG/drE(fp2) = {dGdrE2:.3f}')\n",
    "# print(f'dG/drE(fp3) = {dGdrE3:.3f}')\n",
    "\n",
    "print('\\n')\n",
    "\n",
    "pars = default_pars(wEE=6.4, wEI=4.8, wIE=6.0, wII=1.2, I_ext_E=0.8)\n",
    "x_fp_lc = my_fp(pars, 0.8, 0.8)\n",
    "\n",
    "# dGdrE_lc = get_dGdE(x_fp_lc, **pars)\n",
    "\n",
    "print('For the limit cycle case:')\n",
    "# print(f'dG/drE(fp_lc) = {dGdrE_lc:.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "**SAMPLE OUTPUT**\n",
    "```\n",
    "For the default case:\n",
    "dG/drE(fp1) = -0.650\n",
    "dG/drE(fp2) = 1.519\n",
    "dG/drE(fp3) = -0.706\n",
    "\n",
    "\n",
    "For the limit cycle case:\n",
    "dG/drE(fp_lc) = 0.837\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:21.460869Z",
     "iopub.status.busy": "2021-05-25T01:14:21.459668Z",
     "iopub.status.idle": "2021-05-25T01:14:21.462539Z",
     "shell.execute_reply": "2021-05-25T01:14:21.462103Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "def get_dGdE(fp, tau_E, a_E, theta_E, wEE, wEI, I_ext_E, **other_pars):\n",
    "  \"\"\"\n",
    "  Simulate the Wilson-Cowan equations\n",
    "\n",
    "  Args:\n",
    "    fp   : fixed point (E, I), array\n",
    "    Other arguments are parameters of the Wilson-Cowan model\n",
    "\n",
    "  Returns:\n",
    "    J    : the 2x2 Jacobian matrix\n",
    "  \"\"\"\n",
    "  rE, rI = fp\n",
    "\n",
    "  # Calculate the J[0,0]\n",
    "  dGdrE = (-1 + wEE * dF(wEE * rE - wEI * rI + I_ext_E, a_E, theta_E)) / tau_E\n",
    "\n",
    "  return dGdrE\n",
    "\n",
    "\n",
    "pars = default_pars()\n",
    "x_fp_1 = my_fp(pars, 0.1, 0.1)\n",
    "x_fp_2 = my_fp(pars, 0.3, 0.3)\n",
    "x_fp_3 = my_fp(pars, 0.8, 0.6)\n",
    "\n",
    "dGdrE1 = get_dGdE(x_fp_1, **pars)\n",
    "dGdrE2 = get_dGdE(x_fp_2, **pars)\n",
    "dGdrE3 = get_dGdE(x_fp_3, **pars)\n",
    "\n",
    "print(f'For the default case:')\n",
    "print(f'dG/drE(fp1) = {dGdrE1:.3f}')\n",
    "print(f'dG/drE(fp2) = {dGdrE2:.3f}')\n",
    "print(f'dG/drE(fp3) = {dGdrE3:.3f}')\n",
    "\n",
    "print('\\n')\n",
    "\n",
    "pars = default_pars(wEE=6.4, wEI=4.8, wIE=6.0, wII=1.2, I_ext_E=0.8)\n",
    "x_fp_lc = my_fp(pars, 0.8, 0.8)\n",
    "\n",
    "dGdrE_lc = get_dGdE(x_fp_lc, **pars)\n",
    "\n",
    "print('For the limit cycle case:')\n",
    "print(f'dG/drE(fp_lc) = {dGdrE_lc:.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Nullcline analysis of the ISN\n",
    "\n",
    "Recall that the E nullcline follows\n",
    "\n",
    "\\\\\n",
    "\n",
    "\\begin{align}\n",
    "r_E = F_E(w_{EE}r_E -w_{EI}r_I + I^{\\text{ext}}_E;a_E,\\theta_E). \n",
    "\\end{align}\n",
    "\n",
    "\\\\\n",
    "\n",
    "\n",
    "That is, the firing rate $r_E$ can be a function of $r_I$. Let's take the derivative of $r_E$ over $r_I$, and obtain\n",
    "\n",
    "\\\\\n",
    "\n",
    "\\begin{align}\n",
    "&\\frac{dr_E}{dr_I} = F_E' \\cdot (w_{EE}\\frac{dr_E}{dr_I} -w_{EI}) \\iff \\\\\n",
    "&(1-F_E'w_{EE})\\frac{dr_E}{dr_I} = -F_E' w_{EI} \\iff \\\\\n",
    "&\\frac{dr_E}{dr_I} = \\frac{F_E' w_{EI}}{F_E'w_{EE}-1}.\n",
    "\\end{align}\n",
    "\n",
    "\\\\\n",
    "\n",
    "\n",
    "That is, in the phase plane `rI-rE`-plane, we can obtain the slope along the E nullcline as\n",
    "\n",
    "\\\\\n",
    "\n",
    "\n",
    "$$\\frac{dr_I}{dr_E} = \\frac{F_E'w_{EE}-1}{F_E' w_{EI}} \\qquad (12)$$\n",
    "\n",
    "Similarly, we can obtain the slope along the I nullcline as \n",
    "\n",
    "\\\\\n",
    "\n",
    "$$\\frac{dr_I}{dr_E} = \\frac{F_I'w_{IE}}{F_I' w_{II}+1} \\qquad (13)$$\n",
    "\n",
    "\\\\\n",
    "\n",
    "\n",
    "Then, we can find that $\\Big{(} \\displaystyle{\\frac{dr_I}{dr_E}} \\Big{)}_{\\rm I-nullcline} >0$ in Equation (13).\n",
    "\n",
    "\\\\\n",
    "\n",
    "However, in Equation (12), the sign of $\\Big{(} \\displaystyle{\\frac{dr_I}{dr_E}} \\Big{)}_{\\rm E-nullcline}$ depends on the sign of $(F_E'w_{EE}-1)$. Note that, $(F_E'w_{EE}-1)$ is the same as what we show above (Equation (8)). Therefore, we can have the following results:\n",
    "\n",
    "- $\\Big{(} \\displaystyle{\\frac{dr_I}{dr_E}} \\Big{)}_{\\rm E-nullcline}<0$, **noninhibition-stabilized\n",
    "network (non-ISN) regime**\n",
    "\n",
    "- $\\Big{(} \\displaystyle{\\frac{dr_I}{dr_E}} \\Big{)}_{\\rm E-nullcline}>0$, **inhibition-stabilized\n",
    "network (ISN) regime**\n",
    "\n",
    "\\\\\n",
    "\n",
    "In addition, it is important to point out the following two conclusions: \\\\\n",
    "\n",
    "\n",
    "**Conclusion 1:** The stability of a fixed point can determine the relationship between the slopes Equations (12) and (13). As discussed above, the fixed point is stable when the Jacobian matrix ($J$ in Equation (7)) has two eigenvalues with a negative real part, which indicates a positive determinant of $J$, i.e., $\\text{det}(J)>0$.\n",
    "\n",
    "From the Jacobian matrix definition and from Equations (8-11), we can obtain:\n",
    "\n",
    "$   J=\n",
    "  \\left[ {\\begin{array}{cc}\n",
    "   \\displaystyle{\\frac{1}{\\tau_E}(w_{EE}F_E'-1)} & \\displaystyle{-\\frac{1}{\\tau_E}w_{EI}F_E'}\\\\[1mm]\n",
    "   \\displaystyle {\\frac{1}{\\tau_I}w_{IE}F_I'}& \\displaystyle {\\frac{1}{\\tau_I}(-w_{II}F_I'-1)} \\\\\n",
    "  \\end{array} } \\right] $\n",
    "\n",
    "\\\\\n",
    "\n",
    "Note that, if we let \n",
    "\n",
    "\\\\\n",
    "\n",
    "$   T=\n",
    "  \\left[ {\\begin{array}{cc}\n",
    "   \\displaystyle{\\tau_E} & \\displaystyle{0}\\\\[1mm]\n",
    "   \\displaystyle 0& \\displaystyle \\tau_I \\\\\n",
    "  \\end{array} } \\right] $, \n",
    "$   F=\n",
    "  \\left[ {\\begin{array}{cc}\n",
    "   \\displaystyle{F_E'} & \\displaystyle{0}\\\\[1mm]\n",
    "   \\displaystyle 0& \\displaystyle F_I' \\\\\n",
    "  \\end{array} } \\right] $, and\n",
    "  $   W=\n",
    "  \\left[ {\\begin{array}{cc}\n",
    "   \\displaystyle{w_{EE}} & \\displaystyle{-w_{EI}}\\\\[1mm]\n",
    "   \\displaystyle w_{IE}& \\displaystyle -w_{II} \\\\\n",
    "  \\end{array} } \\right] $\n",
    "\n",
    "\\\\\n",
    "\n",
    "then, using matrix notation, $J=T^{-1}(F W - I)$ where $I$ is the identity matrix, i.e., $I = \\begin{bmatrix} \n",
    "1 & 0 \\\\\n",
    "0 & 1 \n",
    "\\end{bmatrix}.$\n",
    "  \n",
    "\\\\\n",
    "\n",
    "Therefore, $\\det{(J)}=\\det{(T^{-1}(F W - I))}=(\\det{(T^{-1})})(\\det{(F W - I)}).$\n",
    "\n",
    "Since $\\det{(T^{-1})}>0$, as time constants are positive by definition, the sign of $\\det{(J)}$ is the same as the sign of $\\det{(F W - I)}$, and so\n",
    "\n",
    "$$\\det{(FW - I)} = (F_E' w_{EI})(F_I'w_{IE}) - (F_I' w_{II} + 1)(F_E'w_{EE} - 1) > 0.$$\n",
    "\n",
    "\\\\\n",
    "\n",
    "Then, combining this with Equations (12) and (13), we can obtain\n",
    "$$\\frac{\\Big{(} \\displaystyle{\\frac{dr_I}{dr_E}} \\Big{)}_{\\rm I-nullcline}}{\\Big{(} \\displaystyle{\\frac{dr_I}{dr_E}} \\Big{)}_{\\rm E-nullcline}} > 1. $$\n",
    "\n",
    "\n",
    "Therefore, at the stable fixed point, I nullcline has a steeper slope than the E nullcline. \n",
    "\n",
    "\n",
    "**Conclusion 2:** Effect of adding input to the inhibitory population.\n",
    "\n",
    "While adding the input $\\delta I^{\\rm ext}_I$ into the inhibitory population, we can find that the E nullcline (Equation (5)) stays the same, while the I nullcline has a pure left shift: the original I nullcline equation,\n",
    "\n",
    "\\\\\n",
    "\n",
    "\\begin{equation}\n",
    "r_I = F_I(w_{IE}r_E-w_{II}r_I + I^{\\text{ext}}_I ; \\alpha_I, \\theta_I)\n",
    "\\end{equation}\n",
    "\n",
    "\\\\\n",
    "\n",
    "remains true if we take $I^{\\text{ext}}_I \\rightarrow I^{\\text{ext}}_I +\\delta I^{\\rm ext}_I$ and $r_E\\rightarrow r_E'=r_E-\\frac{\\delta I^{\\rm ext}_I}{w_{IE}}$ to obtain\n",
    "\n",
    "\\\\\n",
    "\n",
    "\\begin{equation}\n",
    "r_I = F_I(w_{IE}r_E'-w_{II}r_I + I^{\\text{ext}}_I +\\delta I^{\\rm ext}_I; \\alpha_I, \\theta_I)\n",
    "\\end{equation}\n",
    "\n",
    "\\\\\n",
    "\n",
    "Putting these points together, we obtain the phase plane pictures shown below. After adding input to the inhibitory population, it can be seen in the trajectories above and the phase plane below that, in an **ISN**, $r_I$ will increase first but then decay to the new fixed point in which both $r_I$ and $r_E$ are decreased compared to the original fixed point. However, by adding $\\delta I^{\\rm ext}_I$ into a **non-ISN**, $r_I$ will increase while $r_E$ will decrease."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Interactive Demo: Nullclines of Example **ISN** and **non-ISN**\n",
    "\n",
    "In this interactive widget, we inject excitatory ($I^{\\text{ext}}_I>0$) or inhibitory ($I^{\\text{ext}}_I<0$) drive into the inhibitory population when the system is at its equilibrium (with parameters $w_{EE}=6.4$, $w_{EI}=4.8$, $w_{IE}=6.$, $w_{II}=1.2$, $I_{E}^{\\text{ext}}=0.8$, $\\tau_I = 0.8$, and $I^{\\text{ext}}_I=0$). How does the firing rate of the $I$ population changes with excitatory vs inhibitory drive into the inhibitory population?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:21.490507Z",
     "iopub.status.busy": "2021-05-25T01:14:21.488871Z",
     "iopub.status.idle": "2021-05-25T01:14:21.831091Z",
     "shell.execute_reply": "2021-05-25T01:14:21.830615Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "\n",
    "# @markdown Make sure you execute this cell to enable the widget!\n",
    "\n",
    "pars = default_pars(T=50., dt=0.1)\n",
    "pars['wEE'], pars['wEI'] = 6.4, 4.8\n",
    "pars['wIE'], pars['wII'] = 6.0, 1.2\n",
    "pars['I_ext_E'] = 0.8\n",
    "pars['tau_I'] = 0.8\n",
    "\n",
    "\n",
    "def ISN_I_perturb(dI=0.1):\n",
    "  Lt = len(pars['range_t'])\n",
    "  pars['I_ext_I'] = np.zeros(Lt)\n",
    "  pars['I_ext_I'][int(Lt / 2):] = dI\n",
    "\n",
    "  pars['rE_init'], pars['rI_init'] = 0.6, 0.26\n",
    "  rE, rI = simulate_wc(**pars)\n",
    "\n",
    "  plt.figure(figsize=(8, 1.5))\n",
    "\n",
    "  plt.plot(pars['range_t'], pars['I_ext_I'], 'k')\n",
    "  plt.xlabel('t (ms)')\n",
    "  plt.ylabel(r'$I_I^{\\mathrm{ext}}$')\n",
    "  plt.ylim(pars['I_ext_I'].min() - 0.01, pars['I_ext_I'].max() + 0.01)\n",
    "  plt.show()\n",
    "\n",
    "  plt.figure(figsize=(8, 4.5))\n",
    "  plt.plot(pars['range_t'], rE, 'b', label=r'$r_E$')\n",
    "  plt.plot(pars['range_t'], rE[int(Lt / 2) - 1] * np.ones(Lt), 'b--')\n",
    "  plt.plot(pars['range_t'], rI, 'r', label=r'$r_I$')\n",
    "  plt.plot(pars['range_t'], rI[int(Lt / 2) - 1] * np.ones(Lt), 'r--')\n",
    "  plt.ylim(0, 0.8)\n",
    "  plt.xlabel('t (ms)')\n",
    "  plt.ylabel('Activity')\n",
    "  plt.legend(loc='best')\n",
    "  plt.show()\n",
    "\n",
    "\n",
    "_ = widgets.interact(ISN_I_perturb, dI=(-0.2, 0.21, .05))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:21.842030Z",
     "iopub.status.busy": "2021-05-25T01:14:21.841500Z",
     "iopub.status.idle": "2021-05-25T01:14:21.846454Z",
     "shell.execute_reply": "2021-05-25T01:14:21.845992Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove explanation\n",
    "\n",
    "\"\"\"\n",
    "Discussion:\n",
    "\n",
    "Here we observe a paradoxical effect; if we inject excitatory current to the I\n",
    "population, the r_I goes down, whereas when we inject inhibitory current, the r_I\n",
    "increases. Recall that we inject a constant excitatory current to the E population,\n",
    "which also drives, indirectly, the I population. When Iext>0, the r_I increases\n",
    "but this drives E to a low state, which in turn leads to rI decrease. Whereas,\n",
    "when Iext<0, the effect is negative on I population for a short amount of time,\n",
    "which is sufficient to drive the E population to a high steady state, and then due\n",
    "to E to I connections, the I population activity is increased.\n",
    "\"\"\";"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Bonus 3: Fixed point and working memory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "The input into the neurons measured in the experiment is often very noisy ([links](http://www.scholarpedia.org/article/Stochastic_dynamical_systems)). Here, the noisy synaptic input current is modeled as an Ornstein-Uhlenbeck (OU)process, which has been discussed several times in the previous tutorials.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:21.870884Z",
     "iopub.status.busy": "2021-05-25T01:14:21.868948Z",
     "iopub.status.idle": "2021-05-25T01:14:22.061916Z",
     "shell.execute_reply": "2021-05-25T01:14:22.062351Z"
    }
   },
   "outputs": [],
   "source": [
    "# @markdown Make sure you execute this cell to enable the function my_OU and plot the input current!\n",
    "\n",
    "\n",
    "def my_OU(pars, sig, myseed=False):\n",
    "  \"\"\"\n",
    "  Expects:\n",
    "  pars       : parameter dictionary\n",
    "  sig        : noise amplitute\n",
    "  myseed     : random seed. int or boolean\n",
    "\n",
    "  Returns:\n",
    "  I          : Ornstein-Uhlenbeck input current\n",
    "  \"\"\"\n",
    "\n",
    "  # Retrieve simulation parameters\n",
    "  dt, range_t = pars['dt'], pars['range_t']\n",
    "  Lt = range_t.size\n",
    "  tau_ou = pars['tau_ou']  # [ms]\n",
    "\n",
    "  # set random seed\n",
    "  if myseed:\n",
    "      np.random.seed(seed=myseed)\n",
    "  else:\n",
    "      np.random.seed()\n",
    "\n",
    "  # Initialize\n",
    "  noise = np.random.randn(Lt)\n",
    "  I_ou = np.zeros(Lt)\n",
    "  I_ou[0] = noise[0] * sig\n",
    "\n",
    "  # generate OU\n",
    "  for it in range(Lt-1):\n",
    "      I_ou[it+1] = (I_ou[it]\n",
    "                    + dt / tau_ou * (0. - I_ou[it])\n",
    "                    + np.sqrt(2 * dt / tau_ou) * sig * noise[it + 1])\n",
    "  return I_ou\n",
    "\n",
    "\n",
    "pars = default_pars(T=50)\n",
    "pars['tau_ou'] = 1.  # [ms]\n",
    "sig_ou = 0.1\n",
    "I_ou = my_OU(pars, sig=sig_ou, myseed=2020)\n",
    "plt.figure(figsize=(8, 5.5))\n",
    "plt.plot(pars['range_t'], I_ou, 'b')\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel(r'$I_{\\mathrm{OU}}$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "\n",
    "\n",
    "With the default parameters, the system fluctuates around a resting state with the noisy input.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:22.070427Z",
     "iopub.status.busy": "2021-05-25T01:14:22.069910Z",
     "iopub.status.idle": "2021-05-25T01:14:22.290300Z",
     "shell.execute_reply": "2021-05-25T01:14:22.290737Z"
    }
   },
   "outputs": [],
   "source": [
    "# @markdown Execute this cell to plot activity with noisy input current\n",
    "pars = default_pars(T=100)\n",
    "pars['tau_ou'] = 1.  # [ms]\n",
    "sig_ou = 0.1\n",
    "pars['I_ext_E'] = my_OU(pars, sig=sig_ou, myseed=20201)\n",
    "pars['I_ext_I'] = my_OU(pars, sig=sig_ou, myseed=20202)\n",
    "\n",
    "pars['rE_init'], pars['rI_init'] = 0.1, 0.1\n",
    "rE, rI = simulate_wc(**pars)\n",
    "\n",
    "plt.figure(figsize=(8, 5.5))\n",
    "ax = plt.subplot(111)\n",
    "ax.plot(pars['range_t'], rE, 'b', label='E population')\n",
    "ax.plot(pars['range_t'], rI, 'r', label='I population')\n",
    "ax.set_xlabel('t (ms)')\n",
    "ax.set_ylabel('Activity')\n",
    "ax.legend(loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Interactive Demo: Short pulse induced persistent activity\n",
    "Then, let's use a brief 10-ms positive current to the E population when the system is at its equilibrium. When this amplitude (SE below) is sufficiently large, a persistent activity is produced that outlasts the transient input. What is the firing rate of the persistent activity, and what is the critical input strength? Try to understand the phenomena from the above phase-plane analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:22.300597Z",
     "iopub.status.busy": "2021-05-25T01:14:22.300029Z",
     "iopub.status.idle": "2021-05-25T01:14:22.548172Z",
     "shell.execute_reply": "2021-05-25T01:14:22.547717Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "\n",
    "# @markdown Make sure you execute this cell to enable the widget!\n",
    "\n",
    "\n",
    "def my_inject(pars, t_start, t_lag=10.):\n",
    "  \"\"\"\n",
    "  Expects:\n",
    "  pars       : parameter dictionary\n",
    "  t_start    : pulse starts [ms]\n",
    "  t_lag      : pulse lasts  [ms]\n",
    "\n",
    "  Returns:\n",
    "  I          : extra pulse time\n",
    "  \"\"\"\n",
    "\n",
    "  # Retrieve simulation parameters\n",
    "  dt, range_t = pars['dt'], pars['range_t']\n",
    "  Lt = range_t.size\n",
    "\n",
    "  # Initialize\n",
    "  I = np.zeros(Lt)\n",
    "\n",
    "  # pulse timing\n",
    "  N_start = int(t_start / dt)\n",
    "  N_lag = int(t_lag / dt)\n",
    "  I[N_start:N_start + N_lag] = 1.\n",
    "\n",
    "  return I\n",
    "\n",
    "\n",
    "pars = default_pars(T=100)\n",
    "pars['tau_ou'] = 1.  # [ms]\n",
    "sig_ou = 0.1\n",
    "pars['I_ext_I'] = my_OU(pars, sig=sig_ou, myseed=2021)\n",
    "pars['rE_init'], pars['rI_init'] = 0.1, 0.1\n",
    "\n",
    "# pulse\n",
    "I_pulse = my_inject(pars, t_start=20., t_lag=10.)\n",
    "L_pulse = sum(I_pulse > 0.)\n",
    "\n",
    "\n",
    "def WC_with_pulse(SE=0.):\n",
    "  pars['I_ext_E'] = my_OU(pars, sig=sig_ou, myseed=2022)\n",
    "  pars['I_ext_E'] += SE * I_pulse\n",
    "\n",
    "  rE, rI = simulate_wc(**pars)\n",
    "\n",
    "  plt.figure(figsize=(8, 5.5))\n",
    "  ax = plt.subplot(111)\n",
    "  ax.plot(pars['range_t'], rE, 'b', label='E population')\n",
    "  ax.plot(pars['range_t'], rI, 'r', label='I population')\n",
    "\n",
    "  ax.plot(pars['range_t'][I_pulse > 0.], 1.0*np.ones(L_pulse), 'r', lw=3.)\n",
    "  ax.text(25, 1.05, 'stimulus on', horizontalalignment='center',\n",
    "          verticalalignment='bottom')\n",
    "  ax.set_ylim(-0.03, 1.2)\n",
    "  ax.set_xlabel('t (ms)')\n",
    "  ax.set_ylabel('Activity')\n",
    "  ax.legend(loc='best')\n",
    "  plt.show()\n",
    "\n",
    "\n",
    "_ = widgets.interact(WC_with_pulse, SE=(0.0, 1.0, .05))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:22.555917Z",
     "iopub.status.busy": "2021-05-25T01:14:22.554596Z",
     "iopub.status.idle": "2021-05-25T01:14:22.557818Z",
     "shell.execute_reply": "2021-05-25T01:14:22.557346Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove explanation\n",
    "\n",
    "\"\"\"\n",
    "Discussion:\n",
    "\n",
    "When a system has more than one fixed points, depending on the input strength,\n",
    "the network will settle in one of the fixed points. In this case, we have two\n",
    "fixed points, one of the fixed points corresponds to high activity. So when input\n",
    "drives the network to the high activity fixed points, the network activity will\n",
    "remain there -- it is a stable fixed point. Because the network retains its\n",
    "activity (persistent activity) even after the input has been removed, we can\n",
    "take the persistent activity as working memory.\n",
    "\"\"\";"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "Explore what happened when a second, brief current is applied to the inhibitory population. "
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "W2D4_Tutorial2",
   "provenance": [],
   "toc_visible": true
  },
  "kernel": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
